1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20 import junit.framework.Test;
21 import junit.framework.TestCase;
22 import junit.framework.TestSuite;
23
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;
39
40 public class HighamHall54IntegratorTest
41 extends TestCase {
42
43 public HighamHall54IntegratorTest(String name) {
44 super(name);
45 }
46
47 public void testWrongDerivative() {
48 try {
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 {
56 if (t < -0.5) {
57 throw new DerivativeException("{0}", "oops");
58 } else {
59 throw new DerivativeException(new RuntimeException("oops"));
60 }
61 }
62 public int getDimension() {
63 return 1;
64 }
65 };
66
67 try {
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) {
71
72 }
73
74 try {
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) {
78
79 }
80
81 } catch (Exception e) {
82 fail("wrong exception caught: " + e.getMessage());
83 }
84 }
85
86 public void testMinStep() {
87
88 try {
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 };
94
95 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
96 vecAbsoluteTolerance,
97 vecRelativeTolerance);
98 TestProblemHandler handler = new TestProblemHandler(pb, integ);
99 integ.addStepHandler(handler);
100 integ.integrate(pb,
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) {
107 }
108
109 }
110
111 public void testIncreasingTolerance()
112 throws DerivativeException, IntegratorException {
113
114 int previousCalls = Integer.MAX_VALUE;
115 for (int i = -12; i < -2; ++i) {
116 TestProblem1 pb = new TestProblem1();
117 double minStep = 0;
118 double maxStep = pb.getFinalTime() - pb.getInitialTime();
119 double scalAbsoluteTolerance = Math.pow(10.0, i);
120 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
121
122 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
123 scalAbsoluteTolerance,
124 scalRelativeTolerance);
125 TestProblemHandler handler = new TestProblemHandler(pb, integ);
126 integ.addStepHandler(handler);
127 integ.integrate(pb,
128 pb.getInitialTime(), pb.getInitialState(),
129 pb.getFinalTime(), new double[pb.getDimension()]);
130
131
132
133
134 assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
135 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
136
137 int calls = pb.getCalls();
138 assertEquals(integ.getEvaluations(), calls);
139 assertTrue(calls <= previousCalls);
140 previousCalls = calls;
141
142 }
143
144 }
145
146 public void testBackward()
147 throws DerivativeException, IntegratorException {
148
149 TestProblem5 pb = new TestProblem5();
150 double minStep = 0;
151 double maxStep = pb.getFinalTime() - pb.getInitialTime();
152 double scalAbsoluteTolerance = 1.0e-8;
153 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
154
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()]);
162
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());
167 }
168
169 public void testEvents()
170 throws DerivativeException, IntegratorException {
171
172 TestProblem4 pb = new TestProblem4();
173 double minStep = 0;
174 double maxStep = pb.getFinalTime() - pb.getInitialTime();
175 double scalAbsoluteTolerance = 1.0e-8;
176 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
177
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);
187 }
188 assertEquals(functions.length, integ.getEventHandlers().size());
189 integ.integrate(pb,
190 pb.getInitialTime(), pb.getInitialState(),
191 pb.getFinalTime(), new double[pb.getDimension()]);
192
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());
198
199 }
200
201 public void testEventsErrors() {
202
203 final TestProblem1 pb = new TestProblem1();
204 double minStep = 0;
205 double maxStep = pb.getFinalTime() - pb.getInitialTime();
206 double scalAbsoluteTolerance = 1.0e-8;
207 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
208
209 FirstOrderIntegrator integ =
210 new HighamHall54Integrator(minStep, maxStep,
211 scalAbsoluteTolerance, scalRelativeTolerance);
212 TestProblemHandler handler = new TestProblemHandler(pb, integ);
213 integ.addStepHandler(handler);
214
215 integ.addEventHandler(new EventHandler() {
216 public int eventOccurred(double t, double[] y, boolean increasing) {
217 return EventHandler.CONTINUE;
218 }
219 public double g(double t, double[] y) throws EventException {
220 double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
221 double offset = t - middle;
222 if (offset > 0) {
223 throw new EventException("Evaluation failed for argument = {0}", t);
224 }
225 return offset;
226 }
227 public void resetState(double t, double[] y) {
228 }
229 private static final long serialVersionUID = 935652725339916361L;
230 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
231
232 try {
233 integ.integrate(pb,
234 pb.getInitialTime(), pb.getInitialState(),
235 pb.getFinalTime(), new double[pb.getDimension()]);
236 fail("an exception should have been thrown");
237 } catch (IntegratorException ie) {
238
239 } catch (Exception e) {
240 fail("wrong exception type caught");
241 }
242
243 }
244
245 public void testEventsNoConvergence() {
246
247 final TestProblem1 pb = new TestProblem1();
248 double minStep = 0;
249 double maxStep = pb.getFinalTime() - pb.getInitialTime();
250 double scalAbsoluteTolerance = 1.0e-8;
251 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
252
253 FirstOrderIntegrator integ =
254 new HighamHall54Integrator(minStep, maxStep,
255 scalAbsoluteTolerance, scalRelativeTolerance);
256 TestProblemHandler handler = new TestProblemHandler(pb, integ);
257 integ.addStepHandler(handler);
258
259 integ.addEventHandler(new EventHandler() {
260 public int eventOccurred(double t, double[] y, boolean increasing) {
261 return EventHandler.CONTINUE;
262 }
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);
267 }
268 public void resetState(double t, double[] y) {
269 }
270 private static final long serialVersionUID = 935652725339916361L;
271 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
272
273 try {
274 integ.integrate(pb,
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");
283 }
284
285 }
286
287 public void testSanityChecks() {
288 try {
289 final TestProblem3 pb = new TestProblem3(0.9);
290 double minStep = 0;
291 double maxStep = pb.getFinalTime() - pb.getInitialTime();
292
293 try {
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) {
300
301 }
302
303 try {
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) {
310
311 }
312
313 try {
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) {
320
321 }
322
323 try {
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) {
330
331 }
332
333 try {
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) {
340
341 }
342
343 } catch (Exception e) {
344 fail("wrong exception caught: " + e.getMessage());
345 }
346 }
347
348 public void testKepler()
349 throws DerivativeException, IntegratorException {
350
351 final TestProblem3 pb = new TestProblem3(0.9);
352 double minStep = 0;
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 };
356
357 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
358 vecAbsoluteTolerance,
359 vecRelativeTolerance);
360 integ.addStepHandler(new KeplerHandler(pb));
361 integ.integrate(pb,
362 pb.getInitialTime(), pb.getInitialState(),
363 pb.getFinalTime(), new double[pb.getDimension()]);
364 assertEquals("Higham-Hall 5(4)", integ.getName());
365 }
366
367 private static class KeplerHandler implements StepHandler {
368 public KeplerHandler(TestProblem3 pb) {
369 this.pb = pb;
370 nbSteps = 0;
371 maxError = 0;
372 }
373 public boolean requiresDenseOutput() {
374 return false;
375 }
376 public void reset() {
377 nbSteps = 0;
378 maxError = 0;
379 }
380 public void handleStep(StepInterpolator interpolator,
381 boolean isLast) throws DerivativeException {
382
383 ++nbSteps;
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) {
390 maxError = error;
391 }
392 if (isLast) {
393 assertTrue(maxError < 4e-11);
394 assertTrue(nbSteps < 670);
395 }
396 }
397 private TestProblem3 pb;
398 private int nbSteps;
399 private double maxError;
400 }
401
402 public static Test suite() {
403 return new TestSuite(HighamHall54IntegratorTest.class);
404 }
405
406 }