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