1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
21 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22 import org.apache.commons.math4.legacy.exception.NoBracketingException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24 import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
25 import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
26 import org.apache.commons.math4.legacy.ode.TestProblem1;
27 import org.apache.commons.math4.legacy.ode.TestProblem3;
28 import org.apache.commons.math4.legacy.ode.TestProblem4;
29 import org.apache.commons.math4.legacy.ode.TestProblem5;
30 import org.apache.commons.math4.legacy.ode.TestProblemHandler;
31 import org.apache.commons.math4.legacy.ode.events.EventHandler;
32 import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
33 import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
34 import org.apache.commons.math4.core.jdkmath.JdkMath;
35 import org.junit.Assert;
36 import org.junit.Test;
37
38
39 public class DormandPrince853IntegratorTest {
40
41 @Test
42 public void testMissedEndEvent()
43 throws DimensionMismatchException, NumberIsTooSmallException,
44 MaxCountExceededException, NoBracketingException {
45 final double t0 = 1878250320.0000029;
46 final double tEvent = 1878250379.9999986;
47 final double[] k = { 1.0e-4, 1.0e-5, 1.0e-6 };
48 FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() {
49
50 @Override
51 public int getDimension() {
52 return k.length;
53 }
54
55 @Override
56 public void computeDerivatives(double t, double[] y, double[] yDot) {
57 for (int i = 0; i < y.length; ++i) {
58 yDot[i] = k[i] * y[i];
59 }
60 }
61 };
62
63 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.0, 100.0,
64 1.0e-10, 1.0e-10);
65
66 double[] y0 = new double[k.length];
67 for (int i = 0; i < y0.length; ++i) {
68 y0[i] = i + 1;
69 }
70 double[] y = new double[k.length];
71
72 integrator.setInitialStepSize(60.0);
73 double finalT = integrator.integrate(ode, t0, y0, tEvent, y);
74 Assert.assertEquals(tEvent, finalT, 5.0e-6);
75 for (int i = 0; i < y.length; ++i) {
76 Assert.assertEquals(y0[i] * JdkMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9);
77 }
78
79 integrator.setInitialStepSize(60.0);
80 integrator.addEventHandler(new EventHandler() {
81
82 @Override
83 public void init(double t0, double[] y0, double t) {
84 }
85
86 @Override
87 public void resetState(double t, double[] y) {
88 }
89
90 @Override
91 public double g(double t, double[] y) {
92 return t - tEvent;
93 }
94
95 @Override
96 public Action eventOccurred(double t, double[] y, boolean increasing) {
97 Assert.assertEquals(tEvent, t, 5.0e-6);
98 return Action.CONTINUE;
99 }
100 }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
101 finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y);
102 Assert.assertEquals(tEvent + 120, finalT, 5.0e-6);
103 for (int i = 0; i < y.length; ++i) {
104 Assert.assertEquals(y0[i] * JdkMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9);
105 }
106 }
107
108 @Test(expected=DimensionMismatchException.class)
109 public void testDimensionCheck()
110 throws DimensionMismatchException, NumberIsTooSmallException,
111 MaxCountExceededException, NoBracketingException {
112 TestProblem1 pb = new TestProblem1();
113 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.0, 1.0,
114 1.0e-10, 1.0e-10);
115 integrator.integrate(pb,
116 0.0, new double[pb.getDimension()+10],
117 1.0, new double[pb.getDimension()+10]);
118 Assert.fail("an exception should have been thrown");
119 }
120
121 @Test(expected=NumberIsTooSmallException.class)
122 public void testNullIntervalCheck()
123 throws DimensionMismatchException, NumberIsTooSmallException,
124 MaxCountExceededException, NoBracketingException {
125 TestProblem1 pb = new TestProblem1();
126 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.0, 1.0,
127 1.0e-10, 1.0e-10);
128 integrator.integrate(pb,
129 0.0, new double[pb.getDimension()],
130 0.0, new double[pb.getDimension()]);
131 Assert.fail("an exception should have been thrown");
132 }
133
134 @Test(expected=NumberIsTooSmallException.class)
135 public void testMinStep()
136 throws DimensionMismatchException, NumberIsTooSmallException,
137 MaxCountExceededException, NoBracketingException {
138
139 TestProblem1 pb = new TestProblem1();
140 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
141 double maxStep = pb.getFinalTime() - pb.getInitialTime();
142 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
143 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
144
145 FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
146 vecAbsoluteTolerance,
147 vecRelativeTolerance);
148 TestProblemHandler handler = new TestProblemHandler(pb, integ);
149 integ.addStepHandler(handler);
150 integ.integrate(pb,
151 pb.getInitialTime(), pb.getInitialState(),
152 pb.getFinalTime(), new double[pb.getDimension()]);
153 Assert.fail("an exception should have been thrown");
154 }
155
156 @Test
157 public void testIncreasingTolerance()
158 throws DimensionMismatchException, NumberIsTooSmallException,
159 MaxCountExceededException, NoBracketingException {
160
161 int previousCalls = Integer.MAX_VALUE;
162 AdaptiveStepsizeIntegrator integ =
163 new DormandPrince853Integrator(0, Double.POSITIVE_INFINITY,
164 Double.NaN, Double.NaN);
165 for (int i = -12; i < -2; ++i) {
166 TestProblem1 pb = new TestProblem1();
167 double minStep = 0;
168 double maxStep = pb.getFinalTime() - pb.getInitialTime();
169 double scalAbsoluteTolerance = JdkMath.pow(10.0, i);
170 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
171 integ.setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
172
173 TestProblemHandler handler = new TestProblemHandler(pb, integ);
174 integ.addStepHandler(handler);
175 integ.integrate(pb,
176 pb.getInitialTime(), pb.getInitialState(),
177 pb.getFinalTime(), new double[pb.getDimension()]);
178
179
180
181
182 Assert.assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
183 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
184
185 int calls = pb.getCalls();
186 Assert.assertEquals(integ.getEvaluations(), calls);
187 Assert.assertTrue(calls <= previousCalls);
188 previousCalls = calls;
189 }
190 }
191
192 @Test
193 public void testTooLargeFirstStep()
194 throws DimensionMismatchException, NumberIsTooSmallException,
195 MaxCountExceededException, NoBracketingException {
196
197 AdaptiveStepsizeIntegrator integ =
198 new DormandPrince853Integrator(0, Double.POSITIVE_INFINITY, Double.NaN, Double.NaN);
199 final double start = 0.0;
200 final double end = 0.001;
201 FirstOrderDifferentialEquations equations = new FirstOrderDifferentialEquations() {
202
203 @Override
204 public int getDimension() {
205 return 1;
206 }
207
208 @Override
209 public void computeDerivatives(double t, double[] y, double[] yDot) {
210 Assert.assertTrue(t >= JdkMath.nextAfter(start, Double.NEGATIVE_INFINITY));
211 Assert.assertTrue(t <= JdkMath.nextAfter(end, Double.POSITIVE_INFINITY));
212 yDot[0] = -100.0 * y[0];
213 }
214 };
215
216 integ.setStepSizeControl(0, 1.0, 1.0e-6, 1.0e-8);
217 integ.integrate(equations, start, new double[] { 1.0 }, end, new double[1]);
218 }
219
220 @Test
221 public void testBackward()
222 throws DimensionMismatchException, NumberIsTooSmallException,
223 MaxCountExceededException, NoBracketingException {
224
225 TestProblem5 pb = new TestProblem5();
226 double minStep = 0;
227 double maxStep = pb.getFinalTime() - pb.getInitialTime();
228 double scalAbsoluteTolerance = 1.0e-8;
229 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
230
231 FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
232 scalAbsoluteTolerance,
233 scalRelativeTolerance);
234 TestProblemHandler handler = new TestProblemHandler(pb, integ);
235 integ.addStepHandler(handler);
236 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
237 pb.getFinalTime(), new double[pb.getDimension()]);
238
239 Assert.assertTrue(handler.getLastError() < 1.1e-7);
240 Assert.assertTrue(handler.getMaximalValueError() < 1.1e-7);
241 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
242 Assert.assertEquals("Dormand-Prince 8 (5, 3)", integ.getName());
243 }
244
245 @Test
246 public void testEvents()
247 throws DimensionMismatchException, NumberIsTooSmallException,
248 MaxCountExceededException, NoBracketingException {
249
250 TestProblem4 pb = new TestProblem4();
251 double minStep = 0;
252 double maxStep = pb.getFinalTime() - pb.getInitialTime();
253 double scalAbsoluteTolerance = 1.0e-9;
254 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
255
256 FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
257 scalAbsoluteTolerance,
258 scalRelativeTolerance);
259 TestProblemHandler handler = new TestProblemHandler(pb, integ);
260 integ.addStepHandler(handler);
261 EventHandler[] functions = pb.getEventsHandlers();
262 double convergence = 1.0e-8 * maxStep;
263 for (int l = 0; l < functions.length; ++l) {
264 integ.addEventHandler(functions[l], Double.POSITIVE_INFINITY, convergence, 1000);
265 }
266 Assert.assertEquals(functions.length, integ.getEventHandlers().size());
267 integ.integrate(pb,
268 pb.getInitialTime(), pb.getInitialState(),
269 pb.getFinalTime(), new double[pb.getDimension()]);
270
271 Assert.assertEquals(0, handler.getMaximalValueError(), 2.1e-7);
272 Assert.assertEquals(0, handler.getMaximalTimeError(), convergence);
273 Assert.assertEquals(12.0, handler.getLastTime(), convergence);
274 integ.clearEventHandlers();
275 Assert.assertEquals(0, integ.getEventHandlers().size());
276 }
277
278 @Test
279 public void testKepler()
280 throws DimensionMismatchException, NumberIsTooSmallException,
281 MaxCountExceededException, NoBracketingException {
282
283 final TestProblem3 pb = new TestProblem3(0.9);
284 double minStep = 0;
285 double maxStep = pb.getFinalTime() - pb.getInitialTime();
286 double scalAbsoluteTolerance = 1.0e-8;
287 double scalRelativeTolerance = scalAbsoluteTolerance;
288
289 FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
290 scalAbsoluteTolerance,
291 scalRelativeTolerance);
292 integ.addStepHandler(new KeplerHandler(pb));
293 integ.integrate(pb,
294 pb.getInitialTime(), pb.getInitialState(),
295 pb.getFinalTime(), new double[pb.getDimension()]);
296
297 Assert.assertEquals(integ.getEvaluations(), pb.getCalls());
298 Assert.assertTrue(pb.getCalls() < 3300);
299 }
300
301 @Test
302 public void testVariableSteps()
303 throws DimensionMismatchException, NumberIsTooSmallException,
304 MaxCountExceededException, NoBracketingException {
305
306 final TestProblem3 pb = new TestProblem3(0.9);
307 double minStep = 0;
308 double maxStep = pb.getFinalTime() - pb.getInitialTime();
309 double scalAbsoluteTolerance = 1.0e-8;
310 double scalRelativeTolerance = scalAbsoluteTolerance;
311
312 FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
313 scalAbsoluteTolerance,
314 scalRelativeTolerance);
315 integ.addStepHandler(new VariableHandler());
316 double stopTime = integ.integrate(pb,
317 pb.getInitialTime(), pb.getInitialState(),
318 pb.getFinalTime(), new double[pb.getDimension()]);
319 Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
320 Assert.assertEquals("Dormand-Prince 8 (5, 3)", integ.getName());
321 }
322
323 @Test
324 public void testUnstableDerivative()
325 throws DimensionMismatchException, NumberIsTooSmallException,
326 MaxCountExceededException, NoBracketingException {
327 final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
328 FirstOrderIntegrator integ =
329 new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
330 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
331 double[] y = { Double.NaN };
332 integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
333 Assert.assertEquals(8.0, y[0], 1.0e-12);
334 }
335
336 @Test
337 public void testEventsScheduling() {
338
339 FirstOrderDifferentialEquations sincos = new FirstOrderDifferentialEquations() {
340
341 @Override
342 public int getDimension() {
343 return 2;
344 }
345
346 @Override
347 public void computeDerivatives(double t, double[] y, double[] yDot) {
348 yDot[0] = y[1];
349 yDot[1] = -y[0];
350 }
351 };
352
353 SchedulingChecker sinChecker = new SchedulingChecker(0);
354 SchedulingChecker cosChecker = new SchedulingChecker(1);
355
356 FirstOrderIntegrator integ =
357 new DormandPrince853Integrator(0.001, 1.0, 1.0e-12, 0.0);
358 integ.addEventHandler(sinChecker, 0.01, 1.0e-7, 100);
359 integ.addStepHandler(sinChecker);
360 integ.addEventHandler(cosChecker, 0.01, 1.0e-7, 100);
361 integ.addStepHandler(cosChecker);
362 double t0 = 0.5;
363 double[] y0 = new double[] { JdkMath.sin(t0), JdkMath.cos(t0) };
364 double t = 10.0;
365 double[] y = new double[2];
366 integ.integrate(sincos, t0, y0, t, y);
367 }
368
369 private static final class SchedulingChecker implements StepHandler, EventHandler {
370
371 private int index;
372 private double tMin;
373
374 SchedulingChecker(int index) {
375 this.index = index;
376 }
377
378 @Override
379 public void init(double t0, double[] y0, double t) {
380 tMin = t0;
381 }
382
383 @Override
384 public void handleStep(StepInterpolator interpolator, boolean isLast) {
385 tMin = interpolator.getCurrentTime();
386 }
387
388 @Override
389 public double g(double t, double[] y) {
390
391
392 Assert.assertTrue(t >= tMin);
393 return y[index];
394 }
395
396 @Override
397 public Action eventOccurred(double t, double[] y, boolean increasing) {
398 return Action.RESET_STATE;
399 }
400
401 @Override
402 public void resetState(double t, double[] y) {
403
404 }
405 }
406
407 private static final class KeplerHandler implements StepHandler {
408 KeplerHandler(TestProblem3 pb) {
409 this.pb = pb;
410 }
411 @Override
412 public void init(double t0, double[] y0, double t) {
413 nbSteps = 0;
414 maxError = 0;
415 }
416 @Override
417 public void handleStep(StepInterpolator interpolator, boolean isLast)
418 throws MaxCountExceededException {
419
420 ++nbSteps;
421 for (int a = 1; a < 10; ++a) {
422
423 double prev = interpolator.getPreviousTime();
424 double curr = interpolator.getCurrentTime();
425 double interp = ((10 - a) * prev + a * curr) / 10;
426 interpolator.setInterpolatedTime(interp);
427
428 double[] interpolatedY = interpolator.getInterpolatedState ();
429 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
430 double dx = interpolatedY[0] - theoreticalY[0];
431 double dy = interpolatedY[1] - theoreticalY[1];
432 double error = dx * dx + dy * dy;
433 if (error > maxError) {
434 maxError = error;
435 }
436 }
437 if (isLast) {
438 Assert.assertTrue(maxError < 2.4e-10);
439 Assert.assertTrue(nbSteps < 150);
440 }
441 }
442 private int nbSteps;
443 private double maxError;
444 private TestProblem3 pb;
445 }
446
447 private static final class VariableHandler implements StepHandler {
448 VariableHandler() {
449 firstTime = true;
450 minStep = 0;
451 maxStep = 0;
452 }
453 @Override
454 public void init(double t0, double[] y0, double t) {
455 firstTime = true;
456 minStep = 0;
457 maxStep = 0;
458 }
459 @Override
460 public void handleStep(StepInterpolator interpolator,
461 boolean isLast) {
462
463 double step = JdkMath.abs(interpolator.getCurrentTime()
464 - interpolator.getPreviousTime());
465 if (firstTime) {
466 minStep = JdkMath.abs(step);
467 maxStep = minStep;
468 firstTime = false;
469 } else {
470 if (step < minStep) {
471 minStep = step;
472 }
473 if (step > maxStep) {
474 maxStep = step;
475 }
476 }
477
478 if (isLast) {
479 Assert.assertTrue(minStep < (1.0 / 100.0));
480 Assert.assertTrue(maxStep > (1.0 / 2.0));
481 }
482 }
483 private boolean firstTime = true;
484 private double minStep = 0;
485 private double maxStep = 0;
486 }
487 }