2
* Licensed to the Apache Software Foundation (ASF) under one or more
3
* contributor license agreements. See the NOTICE file distributed with
4
* this work for additional information regarding copyright ownership.
5
* The ASF licenses this file to You under the Apache License, Version 2.0
6
* (the "License"); you may not use this file except in compliance with
7
* the License. You may obtain a copy of the License at
9
* http://www.apache.org/licenses/LICENSE-2.0
11
* Unless required by applicable law or agreed to in writing, software
12
* distributed under the License is distributed on an "AS IS" BASIS,
13
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14
* See the License for the specific language governing permissions and
15
* limitations under the License.
18
package org.apache.commons.math.ode.nonstiff;
20
import junit.framework.Test;
21
import junit.framework.TestCase;
22
import junit.framework.TestSuite;
24
import org.apache.commons.math.ConvergenceException;
25
import org.apache.commons.math.ode.DerivativeException;
26
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
27
import org.apache.commons.math.ode.FirstOrderIntegrator;
28
import org.apache.commons.math.ode.IntegratorException;
29
import org.apache.commons.math.ode.TestProblem1;
30
import org.apache.commons.math.ode.TestProblem3;
31
import org.apache.commons.math.ode.TestProblem4;
32
import org.apache.commons.math.ode.TestProblem5;
33
import org.apache.commons.math.ode.TestProblemHandler;
34
import org.apache.commons.math.ode.events.EventException;
35
import org.apache.commons.math.ode.events.EventHandler;
36
import org.apache.commons.math.ode.nonstiff.HighamHall54Integrator;
37
import org.apache.commons.math.ode.sampling.StepHandler;
38
import org.apache.commons.math.ode.sampling.StepInterpolator;
40
public class HighamHall54IntegratorTest
43
public HighamHall54IntegratorTest(String name) {
47
public void testWrongDerivative() {
49
HighamHall54Integrator integrator =
50
new HighamHall54Integrator(0.0, 1.0, 1.0e-10, 1.0e-10);
51
FirstOrderDifferentialEquations equations =
52
new FirstOrderDifferentialEquations() {
53
private static final long serialVersionUID = -1157081786301178032L;
54
public void computeDerivatives(double t, double[] y, double[] dot)
55
throws DerivativeException {
57
throw new DerivativeException("{0}", "oops");
59
throw new DerivativeException(new RuntimeException("oops"));
62
public int getDimension() {
68
integrator.integrate(equations, -1.0, new double[1], 0.0, new double[1]);
69
fail("an exception should have been thrown");
70
} catch(DerivativeException de) {
75
integrator.integrate(equations, 0.0, new double[1], 1.0, new double[1]);
76
fail("an exception should have been thrown");
77
} catch(DerivativeException de) {
81
} catch (Exception e) {
82
fail("wrong exception caught: " + e.getMessage());
86
public void testMinStep() {
89
TestProblem1 pb = new TestProblem1();
90
double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
91
double maxStep = pb.getFinalTime() - pb.getInitialTime();
92
double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
93
double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
95
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
97
vecRelativeTolerance);
98
TestProblemHandler handler = new TestProblemHandler(pb, integ);
99
integ.addStepHandler(handler);
101
pb.getInitialTime(), pb.getInitialState(),
102
pb.getFinalTime(), new double[pb.getDimension()]);
103
fail("an exception should have been thrown");
104
} catch(DerivativeException de) {
105
fail("wrong exception caught");
106
} catch(IntegratorException ie) {
111
public void testIncreasingTolerance()
112
throws DerivativeException, IntegratorException {
114
int previousCalls = Integer.MAX_VALUE;
115
for (int i = -12; i < -2; ++i) {
116
TestProblem1 pb = new TestProblem1();
118
double maxStep = pb.getFinalTime() - pb.getInitialTime();
119
double scalAbsoluteTolerance = Math.pow(10.0, i);
120
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
122
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
123
scalAbsoluteTolerance,
124
scalRelativeTolerance);
125
TestProblemHandler handler = new TestProblemHandler(pb, integ);
126
integ.addStepHandler(handler);
128
pb.getInitialTime(), pb.getInitialState(),
129
pb.getFinalTime(), new double[pb.getDimension()]);
131
// the 1.3 factor is only valid for this test
132
// and has been obtained from trial and error
133
// there is no general relation between local and global errors
134
assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
135
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
137
int calls = pb.getCalls();
138
assertEquals(integ.getEvaluations(), calls);
139
assertTrue(calls <= previousCalls);
140
previousCalls = calls;
146
public void testBackward()
147
throws DerivativeException, IntegratorException {
149
TestProblem5 pb = new TestProblem5();
151
double maxStep = pb.getFinalTime() - pb.getInitialTime();
152
double scalAbsoluteTolerance = 1.0e-8;
153
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
155
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
156
scalAbsoluteTolerance,
157
scalRelativeTolerance);
158
TestProblemHandler handler = new TestProblemHandler(pb, integ);
159
integ.addStepHandler(handler);
160
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
161
pb.getFinalTime(), new double[pb.getDimension()]);
163
assertTrue(handler.getLastError() < 5.0e-7);
164
assertTrue(handler.getMaximalValueError() < 5.0e-7);
165
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
166
assertEquals("Higham-Hall 5(4)", integ.getName());
169
public void testEvents()
170
throws DerivativeException, IntegratorException {
172
TestProblem4 pb = new TestProblem4();
174
double maxStep = pb.getFinalTime() - pb.getInitialTime();
175
double scalAbsoluteTolerance = 1.0e-8;
176
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
178
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
179
scalAbsoluteTolerance,
180
scalRelativeTolerance);
181
TestProblemHandler handler = new TestProblemHandler(pb, integ);
182
integ.addStepHandler(handler);
183
EventHandler[] functions = pb.getEventsHandlers();
184
for (int l = 0; l < functions.length; ++l) {
185
integ.addEventHandler(functions[l],
186
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
188
assertEquals(functions.length, integ.getEventHandlers().size());
190
pb.getInitialTime(), pb.getInitialState(),
191
pb.getFinalTime(), new double[pb.getDimension()]);
193
assertTrue(handler.getMaximalValueError() < 1.0e-7);
194
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
195
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
196
integ.clearEventHandlers();
197
assertEquals(0, integ.getEventHandlers().size());
201
public void testEventsErrors() {
203
final TestProblem1 pb = new TestProblem1();
205
double maxStep = pb.getFinalTime() - pb.getInitialTime();
206
double scalAbsoluteTolerance = 1.0e-8;
207
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
209
FirstOrderIntegrator integ =
210
new HighamHall54Integrator(minStep, maxStep,
211
scalAbsoluteTolerance, scalRelativeTolerance);
212
TestProblemHandler handler = new TestProblemHandler(pb, integ);
213
integ.addStepHandler(handler);
215
integ.addEventHandler(new EventHandler() {
216
public int eventOccurred(double t, double[] y, boolean increasing) {
217
return EventHandler.CONTINUE;
219
public double g(double t, double[] y) throws EventException {
220
double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
221
double offset = t - middle;
223
throw new EventException("Evaluation failed for argument = {0}", t);
227
public void resetState(double t, double[] y) {
229
private static final long serialVersionUID = 935652725339916361L;
230
}, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
234
pb.getInitialTime(), pb.getInitialState(),
235
pb.getFinalTime(), new double[pb.getDimension()]);
236
fail("an exception should have been thrown");
237
} catch (IntegratorException ie) {
239
} catch (Exception e) {
240
fail("wrong exception type caught");
245
public void testEventsNoConvergence() {
247
final TestProblem1 pb = new TestProblem1();
249
double maxStep = pb.getFinalTime() - pb.getInitialTime();
250
double scalAbsoluteTolerance = 1.0e-8;
251
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
253
FirstOrderIntegrator integ =
254
new HighamHall54Integrator(minStep, maxStep,
255
scalAbsoluteTolerance, scalRelativeTolerance);
256
TestProblemHandler handler = new TestProblemHandler(pb, integ);
257
integ.addStepHandler(handler);
259
integ.addEventHandler(new EventHandler() {
260
public int eventOccurred(double t, double[] y, boolean increasing) {
261
return EventHandler.CONTINUE;
263
public double g(double t, double[] y) {
264
double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
265
double offset = t - middle;
266
return (offset > 0) ? (offset + 0.5) : (offset - 0.5);
268
public void resetState(double t, double[] y) {
270
private static final long serialVersionUID = 935652725339916361L;
271
}, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
275
pb.getInitialTime(), pb.getInitialState(),
276
pb.getFinalTime(), new double[pb.getDimension()]);
277
fail("an exception should have been thrown");
278
} catch (IntegratorException ie) {
279
assertTrue(ie.getCause() != null);
280
assertTrue(ie.getCause() instanceof ConvergenceException);
281
} catch (Exception e) {
282
fail("wrong exception type caught");
287
public void testSanityChecks() {
289
final TestProblem3 pb = new TestProblem3(0.9);
291
double maxStep = pb.getFinalTime() - pb.getInitialTime();
294
FirstOrderIntegrator integ =
295
new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
296
integ.integrate(pb, pb.getInitialTime(), new double[6],
297
pb.getFinalTime(), new double[pb.getDimension()]);
298
fail("an exception should have been thrown");
299
} catch (IntegratorException ie) {
304
FirstOrderIntegrator integ =
305
new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
306
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
307
pb.getFinalTime(), new double[6]);
308
fail("an exception should have been thrown");
309
} catch (IntegratorException ie) {
314
FirstOrderIntegrator integ =
315
new HighamHall54Integrator(minStep, maxStep, new double[2], new double[4]);
316
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
317
pb.getFinalTime(), new double[pb.getDimension()]);
318
fail("an exception should have been thrown");
319
} catch (IntegratorException ie) {
324
FirstOrderIntegrator integ =
325
new HighamHall54Integrator(minStep, maxStep, new double[4], new double[2]);
326
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
327
pb.getFinalTime(), new double[pb.getDimension()]);
328
fail("an exception should have been thrown");
329
} catch (IntegratorException ie) {
334
FirstOrderIntegrator integ =
335
new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
336
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
337
pb.getInitialTime(), new double[pb.getDimension()]);
338
fail("an exception should have been thrown");
339
} catch (IntegratorException ie) {
343
} catch (Exception e) {
344
fail("wrong exception caught: " + e.getMessage());
348
public void testKepler()
349
throws DerivativeException, IntegratorException {
351
final TestProblem3 pb = new TestProblem3(0.9);
353
double maxStep = pb.getFinalTime() - pb.getInitialTime();
354
double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
355
double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
357
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
358
vecAbsoluteTolerance,
359
vecRelativeTolerance);
360
integ.addStepHandler(new KeplerHandler(pb));
362
pb.getInitialTime(), pb.getInitialState(),
363
pb.getFinalTime(), new double[pb.getDimension()]);
364
assertEquals("Higham-Hall 5(4)", integ.getName());
367
private static class KeplerHandler implements StepHandler {
368
public KeplerHandler(TestProblem3 pb) {
373
public boolean requiresDenseOutput() {
376
public void reset() {
380
public void handleStep(StepInterpolator interpolator,
381
boolean isLast) throws DerivativeException {
384
double[] interpolatedY = interpolator.getInterpolatedState();
385
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
386
double dx = interpolatedY[0] - theoreticalY[0];
387
double dy = interpolatedY[1] - theoreticalY[1];
388
double error = dx * dx + dy * dy;
389
if (error > maxError) {
393
assertTrue(maxError < 4e-11);
394
assertTrue(nbSteps < 670);
397
private TestProblem3 pb;
399
private double maxError;
402
public static Test suite() {
403
return new TestSuite(HighamHall54IntegratorTest.class);