1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.ode.nonstiff;
19
20
21 import java.io.ByteArrayInputStream;
22 import java.io.ByteArrayOutputStream;
23 import java.io.IOException;
24 import java.io.ObjectInputStream;
25 import java.io.ObjectOutputStream;
26 import java.util.Random;
27
28 import org.hamcrest.MatcherAssert;
29 import org.hamcrest.Matchers;
30 import org.hipparchus.analysis.UnivariateFunction;
31 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
32 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
33 import org.hipparchus.exception.LocalizedCoreFormats;
34 import org.hipparchus.exception.MathIllegalArgumentException;
35 import org.hipparchus.exception.MathIllegalStateException;
36 import org.hipparchus.ode.DenseOutputModel;
37 import org.hipparchus.ode.ExpandableODE;
38 import org.hipparchus.ode.LocalizedODEFormats;
39 import org.hipparchus.ode.ODEIntegrator;
40 import org.hipparchus.ode.ODEState;
41 import org.hipparchus.ode.ODEStateAndDerivative;
42 import org.hipparchus.ode.OrdinaryDifferentialEquation;
43 import org.hipparchus.ode.SecondaryODE;
44 import org.hipparchus.ode.TestProblem1;
45 import org.hipparchus.ode.TestProblem2;
46 import org.hipparchus.ode.TestProblem3;
47 import org.hipparchus.ode.TestProblem4;
48 import org.hipparchus.ode.TestProblem5;
49 import org.hipparchus.ode.TestProblem6;
50 import org.hipparchus.ode.TestProblemAbstract;
51 import org.hipparchus.ode.TestProblemHandler;
52 import org.hipparchus.ode.events.Action;
53 import org.hipparchus.ode.events.AdaptableInterval;
54 import org.hipparchus.ode.events.ODEEventDetector;
55 import org.hipparchus.ode.events.ODEEventHandler;
56 import org.hipparchus.ode.sampling.ODEStateInterpolator;
57 import org.hipparchus.ode.sampling.ODEStepHandler;
58 import org.hipparchus.ode.sampling.StepInterpolatorTestUtils;
59 import org.hipparchus.util.FastMath;
60 import org.junit.Assert;
61 import org.junit.Test;
62
63 public abstract class RungeKuttaIntegratorAbstractTest {
64
65 protected abstract RungeKuttaIntegrator createIntegrator(double step);
66
67 @Test
68 public abstract void testMissedEndEvent();
69
70 protected void doTestMissedEndEvent(final double epsilonT,
71 final double epsilonY)
72 throws MathIllegalArgumentException, MathIllegalStateException {
73 final double t0 = 1878250320.0000029;
74 final double tEvent = 1878250379.9999986;
75 final double[] k = new double[] { 1.0e-4, 1.0e-5, 1.0e-6 };
76 OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {
77
78 public int getDimension() {
79 return k.length;
80 }
81
82 public double[] computeDerivatives(double t, double[] y) {
83 double[] yDot = new double[k.length];
84 for (int i = 0; i < y.length; ++i) {
85 yDot[i] = k[i] * y[i];
86 }
87 return yDot;
88 }
89 };
90
91 RungeKuttaIntegrator integrator = createIntegrator(60.0);
92
93 double[] y0 = new double[k.length];
94 for (int i = 0; i < y0.length; ++i) {
95 y0[i] = i;
96 }
97
98 ODEStateAndDerivative result = integrator.integrate(new ExpandableODE(ode),
99 new ODEState(t0, y0),
100 tEvent);
101 Assert.assertEquals(tEvent, result.getTime(), epsilonT);
102 double[] y = result.getPrimaryState();
103 for (int i = 0; i < y.length; ++i) {
104 Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (result.getTime() - t0)), y[i], epsilonY);
105 }
106
107 integrator.addEventDetector(new ODEEventDetector() {
108 public AdaptableInterval getMaxCheckInterval() {
109 return s-> Double.POSITIVE_INFINITY;
110 }
111 public int getMaxIterationCount() {
112 return 100;
113 }
114 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
115 return new BracketingNthOrderBrentSolver(0, 1.0e-20, 0, 5);
116 }
117 public double g(ODEStateAndDerivative state) {
118 return state.getTime() - tEvent;
119 }
120 public ODEEventHandler getHandler() {
121 return (state, detector, increasing) -> {
122 Assert.assertEquals(tEvent, state.getTime(), epsilonT);
123 return Action.CONTINUE;
124 };
125 }
126 });
127 result = integrator.integrate(new ExpandableODE(ode),
128 new ODEState(t0, y0),
129 tEvent + 120);
130 Assert.assertEquals(tEvent + 120, result.getTime(), epsilonT);
131 y = result.getPrimaryState();
132 for (int i = 0; i < y.length; ++i) {
133 Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (result.getTime() - t0)), y[i], epsilonY);
134 }
135
136 }
137
138 @Test
139 public abstract void testSanityChecks();
140
141 protected void doTestSanityChecks() {
142 RungeKuttaIntegrator integrator = createIntegrator(0.01);
143 try {
144 TestProblem1 pb = new TestProblem1();
145 integrator.integrate(new ExpandableODE(pb),
146 new ODEState(0.0, new double[pb.getDimension() + 10]),
147 1.0);
148 Assert.fail("an exception should have been thrown");
149 } catch(MathIllegalArgumentException ie) {
150 Assert.assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, ie.getSpecifier());
151 }
152 try {
153 TestProblem1 pb = new TestProblem1();
154 integrator.integrate(new ExpandableODE(pb),
155 new ODEState(0.0, new double[pb.getDimension()]),
156 0.0);
157 Assert.fail("an exception should have been thrown");
158 } catch(MathIllegalArgumentException ie) {
159 Assert.assertEquals(LocalizedODEFormats.TOO_SMALL_INTEGRATION_INTERVAL, ie.getSpecifier());
160 }
161 }
162
163 @Test
164 public abstract void testDecreasingSteps();
165
166 protected void doTestDecreasingSteps(final double safetyValueFactor,
167 final double safetyTimeFactor,
168 final double epsilonT)
169 throws MathIllegalArgumentException, MathIllegalStateException {
170
171 TestProblemAbstract[] allProblems = new TestProblemAbstract[] {
172 new TestProblem1(), new TestProblem2(), new TestProblem3(),
173 new TestProblem4(), new TestProblem5(), new TestProblem6()
174 };
175 for (TestProblemAbstract pb : allProblems) {
176
177 double previousValueError = Double.NaN;
178 double previousTimeError = Double.NaN;
179 for (int i = 4; i < 10; ++i) {
180
181 double step = FastMath.scalb(pb.getFinalTime() - pb.getInitialState().getTime(), -i);
182
183 RungeKuttaIntegrator integ = createIntegrator(step);
184 TestProblemHandler handler = new TestProblemHandler(pb, integ);
185 integ.addStepHandler(handler);
186 double eventTol = 1.0e-6 * step;
187 ODEEventDetector[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY, eventTol, 1000);
188 for (int l = 0; l < functions.length; ++l) {
189 integ.addEventDetector(functions[l]);
190 }
191 Assert.assertEquals(functions.length, integ.getEventDetectors().size());
192 ODEStateAndDerivative stop = integ.integrate(new ExpandableODE(pb),
193 pb.getInitialState(),
194 pb.getFinalTime());
195 if (functions.length == 0) {
196 Assert.assertEquals(pb.getFinalTime(), stop.getTime(), epsilonT);
197 }
198
199 double error = handler.getMaximalValueError();
200 if (i > 4) {
201 Assert.assertTrue(error < FastMath.abs(previousValueError * safetyValueFactor));
202 }
203 previousValueError = error;
204
205 double timeError = handler.getMaximalTimeError();
206
207 double timeTol = FastMath.max(eventTol, FastMath.abs(previousTimeError * safetyTimeFactor));
208 if (i > 4) {
209 MatcherAssert.assertThat(
210 "Problem=" + pb + ", i=" + i + ", step=" + step,
211 timeError,
212 Matchers.lessThanOrEqualTo(timeTol));
213 }
214 previousTimeError = timeError;
215
216 integ.clearEventDetectors();
217 Assert.assertEquals(0, integ.getEventDetectors().size());
218 }
219
220 }
221
222 }
223
224 @Test
225 public abstract void testSmallStep();
226
227 protected void doTestSmallStep(final double epsilonLast,
228 final double epsilonMaxValue,
229 final double epsilonMaxTime,
230 final String name) {
231
232 TestProblem1 pb = new TestProblem1();
233 double step = 0.001 * (pb.getFinalTime() - pb.getInitialState().getTime());
234
235 RungeKuttaIntegrator integ = createIntegrator(step);
236 TestProblemHandler handler = new TestProblemHandler(pb, integ);
237 integ.addStepHandler(handler);
238 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
239
240 Assert.assertEquals(0, handler.getLastError(), epsilonLast);
241 Assert.assertEquals(0, handler.getMaximalValueError(), epsilonMaxValue);
242 Assert.assertEquals(0, handler.getMaximalTimeError(), epsilonMaxTime);
243 Assert.assertEquals(name, integ.getName());
244
245 }
246
247 @Test
248 public abstract void testBigStep();
249
250 protected void doTestBigStep(final double belowLast,
251 final double belowMaxValue,
252 final double epsilonMaxTime,
253 final String name)
254 throws MathIllegalArgumentException, MathIllegalStateException {
255
256 TestProblem1 pb = new TestProblem1();
257 double step = 0.2 * (pb.getFinalTime() - pb.getInitialState().getTime());
258
259 RungeKuttaIntegrator integ = createIntegrator(step);
260 TestProblemHandler handler = new TestProblemHandler(pb, integ);
261 integ.addStepHandler(handler);
262 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
263
264 Assert.assertTrue(handler.getLastError() > belowLast);
265 Assert.assertTrue(handler.getMaximalValueError() > belowMaxValue);
266 Assert.assertEquals(0, handler.getMaximalTimeError(), epsilonMaxTime);
267 Assert.assertEquals(name, integ.getName());
268
269 }
270
271 @Test
272 public abstract void testBackward();
273
274 protected void doTestBackward(final double epsilonLast,
275 final double epsilonMaxValue,
276 final double epsilonMaxTime,
277 final String name)
278 throws MathIllegalArgumentException, MathIllegalStateException {
279
280 TestProblem5 pb = new TestProblem5();
281 double step = FastMath.abs(0.001 * (pb.getFinalTime() - pb.getInitialState().getTime()));
282
283 RungeKuttaIntegrator integ = createIntegrator(step);
284 TestProblemHandler handler = new TestProblemHandler(pb, integ);
285 integ.addStepHandler(handler);
286 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
287
288 Assert.assertEquals(0, handler.getLastError(), epsilonLast);
289 Assert.assertEquals(0, handler.getMaximalValueError(), epsilonMaxValue);
290 Assert.assertEquals(0, handler.getMaximalTimeError(), epsilonMaxTime);
291 Assert.assertEquals(name, integ.getName());
292
293 }
294
295 @Test
296 public abstract void testKepler();
297
298 protected void doTestKepler(double expectedMaxError, double epsilon)
299 throws MathIllegalArgumentException, MathIllegalStateException {
300
301 final TestProblem3 pb = new TestProblem3(0.9);
302 double step = 0.0003 * (pb.getFinalTime() - pb.getInitialState().getTime());
303
304 RungeKuttaIntegrator integ = createIntegrator(step);
305 integ.addStepHandler(new KeplerHandler(pb, expectedMaxError, epsilon));
306 final ExpandableODE expandable = new ExpandableODE(pb);
307 Assert.assertSame(pb, expandable.getPrimary());
308 integ.integrate(expandable, pb.getInitialState(), pb.getFinalTime());
309 }
310
311 private static class KeplerHandler implements ODEStepHandler {
312 private double maxError;
313 private final TestProblem3 pb;
314 private final double expectedMaxError;
315 private final double epsilon;
316 public KeplerHandler(TestProblem3 pb, double expectedMaxError, double epsilon) {
317 this.pb = pb;
318 this.expectedMaxError = expectedMaxError;
319 this.epsilon = epsilon;
320 this.maxError = 0;
321 }
322 public void init(ODEStateAndDerivative state0, double t) {
323 maxError = 0;
324 }
325 public void handleStep(ODEStateInterpolator interpolator) {
326
327 ODEStateAndDerivative current = interpolator.getCurrentState();
328 double[] theoreticalY = pb.computeTheoreticalState(current.getTime());
329 double dx = current.getPrimaryState()[0] - theoreticalY[0];
330 double dy = current.getPrimaryState()[1] - theoreticalY[1];
331 maxError = FastMath.max(maxError, dx * dx + dy * dy);
332 }
333 public void finish(ODEStateAndDerivative finalState) {
334 Assert.assertEquals(expectedMaxError, maxError, epsilon);
335 }
336 }
337
338 @Test
339 public abstract void testStepSize();
340
341 protected void doTestStepSize(final double epsilon)
342 throws MathIllegalArgumentException, MathIllegalStateException {
343 final double finalTime = 5.0;
344 final double step = 1.23456;
345 RungeKuttaIntegrator integ = createIntegrator(step);
346 integ.addStepHandler(new ODEStepHandler() {
347 public void handleStep(ODEStateInterpolator interpolator) {
348 if (interpolator.getCurrentState().getTime() < finalTime - 0.001) {
349 Assert.assertEquals(step,
350 interpolator.getCurrentState().getTime() - interpolator.getPreviousState().getTime(),
351 epsilon);
352 }
353 }
354 });
355 integ.integrate(new ExpandableODE(new OrdinaryDifferentialEquation() {
356 public double[] computeDerivatives(double t, double[] y) {
357 return new double[] { 1.0 };
358 }
359 public int getDimension() {
360 return 1;
361 }
362 }), new ODEState(0, new double[1]), finalTime);
363 }
364
365 @Test
366 public abstract void testSingleStep();
367
368 protected void doTestSingleStep(final double epsilon) {
369
370 final TestProblem3 pb = new TestProblem3(0.9);
371 double h = 0.0003 * (pb.getFinalTime() - pb.getInitialState().getTime());
372
373 RungeKuttaIntegrator integ = createIntegrator(Double.NaN);
374 double t = pb.getInitialState().getTime();
375 double[] y = pb.getInitialState().getPrimaryState();
376 for (int i = 0; i < 100; ++i) {
377 y = integ.singleStep(pb, t, y, t + h);
378 t += h;
379 }
380 double[] yth = pb.computeTheoreticalState(t);
381 double dx = y[0] - yth[0];
382 double dy = y[1] - yth[1];
383 double error = dx * dx + dy * dy;
384 Assert.assertEquals(0.0, error, epsilon);
385 }
386
387 @Test
388 public abstract void testTooLargeFirstStep();
389
390 protected void doTestTooLargeFirstStep() {
391
392 RungeKuttaIntegrator integ = createIntegrator(0.5);
393 final double t0 = 0;
394 final double[] y0 = new double[] { 1.0 };
395 final double t = 0.001;
396 OrdinaryDifferentialEquation equations = new OrdinaryDifferentialEquation() {
397
398 public int getDimension() {
399 return 1;
400 }
401
402 public double[] computeDerivatives(double t, double[] y) {
403 Assert.assertTrue(t >= FastMath.nextAfter(t0, Double.NEGATIVE_INFINITY));
404 Assert.assertTrue(t <= FastMath.nextAfter(t, Double.POSITIVE_INFINITY));
405 return new double[] { -100 * y[0] };
406 }
407
408 };
409
410 integ.integrate(new ExpandableODE(equations), new ODEState(t0, y0), t);
411
412 }
413
414 @Test
415 public abstract void testUnstableDerivative();
416
417 protected void doTestUnstableDerivative(double epsilon) {
418 final StepProblem stepProblem = new StepProblem(s -> 999.0, 1.0e+12, 1000000, 0.0, 1.0, 2.0).
419 withMaxCheck(1.0).
420 withMaxIter(1000).
421 withThreshold(1.0e-12);
422 Assert.assertEquals(1.0, stepProblem.getMaxCheckInterval().currentInterval(null), 1.0e-15);
423 Assert.assertEquals(1000, stepProblem.getMaxIterationCount());
424 Assert.assertEquals(1.0e-12, stepProblem.getSolver().getAbsoluteAccuracy(), 1.0e-25);
425 Assert.assertNotNull(stepProblem.getHandler());
426 RungeKuttaIntegrator integ = createIntegrator(0.3);
427 integ.addEventDetector(stepProblem);
428 ODEStateAndDerivative result = integ.integrate(new ExpandableODE(stepProblem),
429 new ODEState(0, new double[1]),
430 10.0);
431 Assert.assertEquals(8.0, result.getPrimaryState()[0], epsilon);
432 }
433
434 @Test
435 public abstract void testDerivativesConsistency();
436
437 protected void doTestDerivativesConsistency(double epsilon) {
438 TestProblem3 pb = new TestProblem3();
439 double step = 0.001 * (pb.getFinalTime() - pb.getInitialState().getTime());
440 RungeKuttaIntegrator integ = createIntegrator(step);
441 StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.001, 1.0e-10);
442 }
443
444 @Test
445 public abstract void testSecondaryEquations();
446
447 protected void doTestSecondaryEquations(final double epsilonSinCos,
448 final double epsilonLinear) {
449 OrdinaryDifferentialEquation sinCos = new OrdinaryDifferentialEquation() {
450
451 @Override
452 public int getDimension() {
453 return 2;
454 }
455
456 @Override
457 public double[] computeDerivatives(double t, double[] y) {
458
459
460
461
462
463 return new double[] { 0.5 * y[1], -0.5 * y[0] };
464 }
465
466 };
467
468 SecondaryODE linear = new SecondaryODE() {
469
470 @Override
471 public int getDimension() {
472 return 1;
473 }
474
475 @Override
476 public double[] computeDerivatives(double t, double[] primary, double[] primaryDot, double[] secondary) {
477 for (int i = 0; i < primaryDot.length; ++i) {
478
479
480
481
482 primaryDot[i] *= 2;
483 }
484 return new double[] { -1 };
485 }
486
487 };
488
489 ExpandableODE expandable = new ExpandableODE(sinCos);
490 expandable.addSecondaryEquations(linear);
491
492 ODEIntegrator integrator = createIntegrator(0.001);
493 final double[] max = new double[2];
494 integrator.addStepHandler(new ODEStepHandler() {
495 @Override
496 public void handleStep(ODEStateInterpolator interpolator) {
497 for (int i = 0; i <= 10; ++i) {
498 double tPrev = interpolator.getPreviousState().getTime();
499 double tCurr = interpolator.getCurrentState().getTime();
500 double t = (tPrev * (10 - i) + tCurr * i) / 10;
501 ODEStateAndDerivative state = interpolator.getInterpolatedState(t);
502 Assert.assertEquals(2, state.getPrimaryStateDimension());
503 Assert.assertEquals(1, state.getNumberOfSecondaryStates());
504 Assert.assertEquals(2, state.getSecondaryStateDimension(0));
505 Assert.assertEquals(1, state.getSecondaryStateDimension(1));
506 Assert.assertEquals(3, state.getCompleteStateDimension());
507 max[0] = FastMath.max(max[0],
508 FastMath.abs(FastMath.sin(t) - state.getPrimaryState()[0]));
509 max[0] = FastMath.max(max[0],
510 FastMath.abs(FastMath.cos(t) - state.getPrimaryState()[1]));
511 max[1] = FastMath.max(max[1],
512 FastMath.abs(1 - t - state.getSecondaryState(1)[0]));
513 }
514 }
515 });
516
517 double[] primary0 = new double[] { 0.0, 1.0 };
518 double[][] secondary0 = new double[][] { { 1.0 } };
519 ODEState initialState = new ODEState(0.0, primary0, secondary0);
520
521 ODEStateAndDerivative finalState =
522 integrator.integrate(expandable, initialState, 10.0);
523 Assert.assertEquals(10.0, finalState.getTime(), 1.0e-12);
524 Assert.assertEquals(0, max[0], epsilonSinCos);
525 Assert.assertEquals(0, max[1], epsilonLinear);
526
527 }
528
529 @Test
530 public void testNaNAppearing() {
531 try {
532 ODEIntegrator integ = createIntegrator(0.3);
533 integ.integrate(new OrdinaryDifferentialEquation() {
534 public int getDimension() {
535 return 1;
536 }
537 public double[] computeDerivatives(double t, double[] y) {
538 return new double[] { FastMath.log(t) };
539 }
540 }, new ODEState(1.0, new double[] { 1.0 }), -1.0);
541 Assert.fail("an exception should have been thrown");
542 } catch (MathIllegalStateException mise) {
543 Assert.assertEquals(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION, mise.getSpecifier());
544 Assert.assertTrue(((Double) mise.getParts()[0]).doubleValue() <= 0.0);
545 }
546 }
547
548 @Test
549 public abstract void testSerialization();
550
551 protected void doTestSerialization(int expectedSize, double tolerance) {
552 try {
553 TestProblem3 pb = new TestProblem3(0.9);
554 double h = 0.0003 * (pb.getFinalTime() - pb.getInitialState().getTime());
555 RungeKuttaIntegrator integ = createIntegrator(h);
556
557 integ.addStepHandler(new DenseOutputModel());
558 integ.integrate(pb, pb.getInitialState(), pb.getFinalTime());
559
560 ByteArrayOutputStream bos = new ByteArrayOutputStream();
561 ObjectOutputStream oos = new ObjectOutputStream(bos);
562 for (ODEStepHandler handler : integ.getStepHandlers()) {
563 oos.writeObject(handler);
564 }
565
566 Assert.assertTrue("size = " + bos.size (), bos.size () > 9 * expectedSize / 10);
567 Assert.assertTrue("size = " + bos.size (), bos.size () < 11 * expectedSize / 10);
568
569 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
570 ObjectInputStream ois = new ObjectInputStream(bis);
571 DenseOutputModel cm = (DenseOutputModel) ois.readObject();
572
573 Random random = new Random(347588535632l);
574 double maxError = 0.0;
575 for (int i = 0; i < 1000; ++i) {
576 double r = random.nextDouble();
577 double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
578 double[] interpolatedY = cm.getInterpolatedState(time).getPrimaryState();
579 double[] theoreticalY = pb.computeTheoreticalState(time);
580 double dx = interpolatedY[0] - theoreticalY[0];
581 double dy = interpolatedY[1] - theoreticalY[1];
582 double error = dx * dx + dy * dy;
583 if (error > maxError) {
584 maxError = error;
585 }
586 }
587
588 Assert.assertEquals(0, maxError, tolerance);
589
590 } catch (IOException | ClassNotFoundException e) {
591 Assert.fail(e.getLocalizedMessage());
592 }
593
594 }
595
596 @Test
597 public void testIssue250() {
598 final double defaultStep = 60.;
599 RungeKuttaIntegrator integrator = createIntegrator(defaultStep);
600 Assert.assertEquals(defaultStep, integrator.getDefaultStep(), 0.);
601 }
602
603 }