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.hipparchus.analysis.UnivariateFunction;
21 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
22 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
23 import org.hipparchus.exception.MathIllegalArgumentException;
24 import org.hipparchus.exception.MathIllegalStateException;
25 import org.hipparchus.geometry.euclidean.threed.Rotation;
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.hipparchus.ode.ExpandableODE;
28 import org.hipparchus.ode.LocalizedODEFormats;
29 import org.hipparchus.ode.ODEIntegrator;
30 import org.hipparchus.ode.ODEJacobiansProvider;
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.TestProblem3;
37 import org.hipparchus.ode.TestProblem4;
38 import org.hipparchus.ode.TestProblem5;
39 import org.hipparchus.ode.TestProblem7;
40 import org.hipparchus.ode.TestProblem8;
41 import org.hipparchus.ode.TestProblemHandler;
42 import org.hipparchus.ode.VariationalEquation;
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.events.ODEStepEndHandler;
48 import org.hipparchus.ode.sampling.ODEStateInterpolator;
49 import org.hipparchus.ode.sampling.ODEStepHandler;
50 import org.hipparchus.util.CombinatoricsUtils;
51 import org.hipparchus.util.FastMath;
52 import org.hipparchus.util.SinCos;
53 import org.junit.jupiter.api.Test;
54
55 import java.lang.reflect.InvocationTargetException;
56 import java.lang.reflect.Method;
57 import java.util.ArrayList;
58 import java.util.Arrays;
59 import java.util.List;
60 import java.util.stream.Stream;
61
62 import static org.junit.jupiter.api.Assertions.assertEquals;
63 import static org.junit.jupiter.api.Assertions.assertNotNull;
64 import static org.junit.jupiter.api.Assertions.assertSame;
65 import static org.junit.jupiter.api.Assertions.assertTrue;
66 import static org.junit.jupiter.api.Assertions.fail;
67
68
69 public abstract class EmbeddedRungeKuttaIntegratorAbstractTest {
70
71 protected abstract EmbeddedRungeKuttaIntegrator
72 createIntegrator(final double minStep, final double maxStep,
73 final double scalAbsoluteTolerance, final double scalRelativeTolerance);
74
75 protected abstract EmbeddedRungeKuttaIntegrator
76 createIntegrator(final double minStep, final double maxStep,
77 final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
78
79 @Test
80 public abstract void testForwardBackwardExceptions();
81
82 protected void doTestForwardBackwardExceptions() {
83 OrdinaryDifferentialEquation equations = new OrdinaryDifferentialEquation() {
84
85 public int getDimension() {
86 return 1;
87 }
88
89 public double[] computeDerivatives(double t, double[] y) {
90 if (t < -0.5) {
91 throw new LocalException();
92 } else {
93 throw new RuntimeException("oops");
94 }
95 }
96 };
97
98 EmbeddedRungeKuttaIntegrator integrator = createIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
99
100 try {
101 integrator.integrate(new ExpandableODE(equations),
102 new ODEState(-1, new double[1]),
103 0);
104 fail("an exception should have been thrown");
105 } catch(LocalException de) {
106
107 }
108
109 try {
110 integrator.integrate(new ExpandableODE(equations),
111 new ODEState(0, new double[1]),
112 1);
113 fail("an exception should have been thrown");
114 } catch(RuntimeException de) {
115
116 }
117 }
118
119 protected static class LocalException extends RuntimeException {
120 private static final long serialVersionUID = 20151208L;
121 }
122
123 @Test
124 public void testMinStep() {
125 try {
126 TestProblem1 pb = new TestProblem1();
127 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialState().getTime());
128 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
129 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
130 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
131
132 ODEIntegrator integ = createIntegrator(minStep, maxStep,
133 vecAbsoluteTolerance, vecRelativeTolerance);
134 TestProblemHandler handler = new TestProblemHandler(pb, integ);
135 integ.addStepHandler(handler);
136 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
137 fail("an exception should have been thrown");
138 } catch (MathIllegalArgumentException miae) {
139 assertEquals(LocalizedODEFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
140 miae.getSpecifier());
141 }
142
143 }
144
145 @Test
146 public abstract void testIncreasingTolerance();
147
148 protected void doTestIncreasingTolerance(double factor, double epsilon) {
149
150 int previousCalls = Integer.MAX_VALUE;
151 for (int i = -12; i < -2; ++i) {
152 TestProblem1 pb = new TestProblem1();
153 double minStep = 0;
154 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
155 double scalAbsoluteTolerance = FastMath.pow(10.0, i);
156 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
157
158 ODEIntegrator integ = createIntegrator(minStep, maxStep,
159 scalAbsoluteTolerance, scalRelativeTolerance);
160 TestProblemHandler handler = new TestProblemHandler(pb, integ);
161 integ.addStepHandler(handler);
162 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
163
164 assertTrue(handler.getMaximalValueError() < (factor * scalAbsoluteTolerance));
165 assertEquals(0, handler.getMaximalTimeError(), epsilon);
166
167 int calls = pb.getCalls();
168 assertEquals(integ.getEvaluations(), calls);
169 assertTrue(calls <= previousCalls);
170 previousCalls = calls;
171
172 }
173
174 }
175
176 @Test
177 public abstract void testEvents();
178
179 protected void doTestEvents(final double epsilonMaxValue, final String name) {
180
181 TestProblem4 pb = new TestProblem4();
182 double minStep = 0;
183 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
184 double scalAbsoluteTolerance = 1.0e-8;
185 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
186
187 ODEIntegrator integ = createIntegrator(minStep, maxStep,
188 scalAbsoluteTolerance, scalRelativeTolerance);
189 TestProblemHandler handler = new TestProblemHandler(pb, integ);
190 integ.addStepHandler(handler);
191 double convergence = 1.0e-8 * maxStep;
192 ODEEventDetector[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY, convergence, 1000);
193 for (int l = 0; l < functions.length; ++l) {
194 integ.addEventDetector(functions[l]);
195 }
196 List<ODEEventDetector> detectors = new ArrayList<>(integ.getEventDetectors());
197 assertEquals(functions.length, detectors.size());
198
199 for (int i = 0; i < detectors.size(); ++i) {
200 assertSame(functions[i], detectors.get(i).getHandler());
201 assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null, true), 1.0);
202 assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy(), 1.0e-15 * convergence);
203 assertEquals(1000, detectors.get(i).getMaxIterationCount());
204 }
205
206 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
207
208 assertEquals(0, handler.getMaximalValueError(), epsilonMaxValue);
209 assertEquals(0, handler.getMaximalTimeError(), convergence);
210 assertEquals(12.0, handler.getLastTime(), convergence);
211 assertEquals(name, integ.getName());
212 integ.clearEventDetectors();
213 assertEquals(0, integ.getEventDetectors().size());
214
215 }
216
217 @Test
218 public abstract void testStepEnd();
219
220 protected void doTestStepEnd(final int expectedCount, final String name) {
221 TestProblem4 pb = new TestProblem4();
222 double minStep = 0;
223 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
224 double scalAbsoluteTolerance = 1.0e-8;
225 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
226
227 ODEIntegrator integ = createIntegrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
228 double convergence = 1.0e-8 * maxStep;
229 ODEEventDetector[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY, convergence, 1000);
230 for (int l = 0; l < functions.length; ++l) {
231 integ.addEventDetector(functions[l]);
232 }
233 List<ODEEventDetector> detectors = new ArrayList<>(integ.getEventDetectors());
234 assertEquals(functions.length, detectors.size());
235
236 for (int i = 0; i < detectors.size(); ++i) {
237 assertSame(functions[i], detectors.get(i).getHandler());
238 assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null, true), 1.0);
239 assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy(), 1.0e-15 * convergence);
240 assertEquals(1000, detectors.get(i).getMaxIterationCount());
241 }
242
243 final StepCounter counter = new StepCounter(expectedCount + 10, Action.STOP);
244 integ.addStepEndHandler(counter);
245 assertEquals(1, integ.getStepEndHandlers().size());
246 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
247
248 assertEquals(expectedCount, counter.count);
249 assertEquals(name, integ.getName());
250 integ.clearEventDetectors();
251 assertEquals(0, integ.getEventDetectors().size());
252 integ.clearStepEndHandlers();
253 assertEquals(0, integ.getStepEndHandlers().size());
254 }
255
256 @Test
257 public abstract void testStopAfterStep();
258
259 protected void doTestStopAfterStep(final int count, final double expectedTime) {
260 TestProblem4 pb = new TestProblem4();
261 double minStep = 0;
262 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
263 double scalAbsoluteTolerance = 1.0e-8;
264 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
265
266 ODEIntegrator integ = createIntegrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
267 double convergence = 1.0e-8 * maxStep;
268 ODEEventDetector[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY, convergence, 1000);
269 for (int l = 0; l < functions.length; ++l) {
270 integ.addEventDetector(functions[l]);
271 }
272 List<ODEEventDetector> detectors = new ArrayList<>(integ.getEventDetectors());
273 assertEquals(functions.length, detectors.size());
274
275 for (int i = 0; i < detectors.size(); ++i) {
276 assertSame(functions[i], detectors.get(i).getHandler());
277 assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null, true), 1.0);
278 assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy(), 1.0e-15 * convergence);
279 assertEquals(1000, detectors.get(i).getMaxIterationCount());
280 }
281
282 final StepCounter counter = new StepCounter(count, Action.STOP);
283 integ.addStepEndHandler(counter);
284 assertEquals(1, integ.getStepEndHandlers().size());
285 ODEStateAndDerivative finalState = integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
286
287 assertEquals(count, counter.count);
288 assertEquals(expectedTime, finalState.getTime(), 1.0e-6);
289
290 }
291
292 @Test
293 public abstract void testResetAfterStep();
294
295 protected void doTestResetAfterStep(final int resetCount, final int expectedCount) {
296 TestProblem4 pb = new TestProblem4();
297 double minStep = 0;
298 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
299 double scalAbsoluteTolerance = 1.0e-8;
300 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
301
302 ODEIntegrator integ = createIntegrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
303 double convergence = 1.0e-8 * maxStep;
304 ODEEventDetector[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY, convergence, 1000);
305 for (int l = 0; l < functions.length; ++l) {
306 integ.addEventDetector(functions[l]);
307 }
308 List<ODEEventDetector> detectors = new ArrayList<>(integ.getEventDetectors());
309 assertEquals(functions.length, detectors.size());
310
311 for (int i = 0; i < detectors.size(); ++i) {
312 assertSame(functions[i], detectors.get(i).getHandler());
313 assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null, true), 1.0);
314 assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy(), 1.0e-15 * convergence);
315 assertEquals(1000, detectors.get(i).getMaxIterationCount());
316 }
317
318 final StepCounter counter = new StepCounter(resetCount, Action.RESET_STATE);
319 integ.addStepEndHandler(counter);
320 assertEquals(1, integ.getStepEndHandlers().size());
321 ODEStateAndDerivative finalState = integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
322
323 assertEquals(expectedCount, counter.count);
324 assertEquals(12.0, finalState.getTime(), 1.0e-6);
325 for (int i = 0; i < finalState.getPrimaryStateDimension(); ++i) {
326 assertEquals(0.0, finalState.getPrimaryState()[i], 1.0e-15);
327 assertEquals(0.0, finalState.getPrimaryDerivative()[i], 1.0e-15);
328 }
329
330 }
331
332 private static class StepCounter implements ODEStepEndHandler {
333 final int max;
334 final Action actionAtMax;
335 int count;
336 StepCounter(final int max, final Action actionAtMax) {
337 this.max = max;
338 this.actionAtMax = actionAtMax;
339 this.count = 0;
340 }
341 public Action stepEndOccurred(final ODEStateAndDerivative state, final boolean forward) {
342 return ++count == max ? actionAtMax : Action.CONTINUE;
343 }
344 public ODEState resetState(ODEStateAndDerivative state) {
345 return new ODEState(state.getTime(),
346 new double[state.getPrimaryStateDimension()]);
347 }
348 }
349
350 @Test
351 public void testEventsErrors() {
352 try {
353 final TestProblem1 pb = new TestProblem1();
354 double minStep = 0;
355 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
356 double scalAbsoluteTolerance = 1.0e-8;
357 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
358
359 ODEIntegrator integ = createIntegrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
360 TestProblemHandler handler = new TestProblemHandler(pb, integ);
361 integ.addStepHandler(handler);
362
363 integ.addEventDetector(new ODEEventDetector() {
364 public AdaptableInterval getMaxCheckInterval() {
365 return (s, isForward)-> Double.POSITIVE_INFINITY;
366 }
367 public int getMaxIterationCount() {
368 return 1000;
369 }
370 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
371 return new BracketingNthOrderBrentSolver(0, 1.0e-8 * maxStep, 0, 5);
372 }
373 public ODEEventHandler getHandler() {
374 return (state, detector, increasing) -> Action.CONTINUE;
375 }
376 public double g(ODEStateAndDerivative state) {
377 double middle = 0.5 * (pb.getInitialState().getTime() + pb.getFinalTime());
378 double offset = state.getTime() - middle;
379 if (offset > 0) {
380 throw new LocalException();
381 }
382 return offset;
383 }
384 });
385
386 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
387 } catch (LocalException le) {
388
389 }
390
391 }
392
393 @Test
394 public void testEventsNoConvergence() {
395
396 final TestProblem1 pb = new TestProblem1();
397 double minStep = 0;
398 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
399 double scalAbsoluteTolerance = 1.0e-8;
400 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
401
402 ODEIntegrator integ = createIntegrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
403 TestProblemHandler handler = new TestProblemHandler(pb, integ);
404 integ.addStepHandler(handler);
405
406 integ.addEventDetector(new ODEEventDetector() {
407 public AdaptableInterval getMaxCheckInterval() {
408 return (s, isForward)-> Double.POSITIVE_INFINITY;
409 }
410 public int getMaxIterationCount() {
411 return 3;
412 }
413 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
414 return new BracketingNthOrderBrentSolver(0, 1.0e-8 * maxStep, 0, 5);
415 }
416 public ODEEventHandler getHandler() {
417 return (state, detector, increasing) -> Action.CONTINUE;
418 }
419 public double g(ODEStateAndDerivative state) {
420 double middle = 0.5 * (pb.getInitialState().getTime() + pb.getFinalTime());
421 double offset = state.getTime() - middle;
422 return (offset > 0) ? offset + 0.5 : offset - 0.5;
423 }
424 });
425
426 try {
427 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
428 fail("an exception should have been thrown");
429 } catch (MathIllegalStateException mcee) {
430
431 }
432
433 }
434
435 @Test
436 public void testSanityChecks() {
437 TestProblem3 pb = new TestProblem3();
438 try {
439 EmbeddedRungeKuttaIntegrator integrator = createIntegrator(0,
440 pb.getFinalTime() - pb.getInitialState().getTime(),
441 new double[4], new double[4]);
442 integrator.integrate(new ExpandableODE(pb),
443 new ODEState(pb.getInitialState().getTime(), new double[6]),
444 pb.getFinalTime());
445 fail("an exception should have been thrown");
446 } catch(MathIllegalArgumentException ie) {
447 }
448 try {
449 EmbeddedRungeKuttaIntegrator integrator =
450 createIntegrator(0,
451 pb.getFinalTime() - pb.getInitialState().getTime(),
452 new double[2], new double[4]);
453 integrator.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
454 fail("an exception should have been thrown");
455 } catch(MathIllegalArgumentException ie) {
456 }
457 try {
458 EmbeddedRungeKuttaIntegrator integrator =
459 createIntegrator(0,
460 pb.getFinalTime() - pb.getInitialState().getTime(),
461 new double[4], new double[4]);
462 integrator.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getInitialState().getTime());
463 fail("an exception should have been thrown");
464 } catch(MathIllegalArgumentException ie) {
465 }
466 }
467
468 @Test
469 public void testNullIntervalCheck()
470 throws MathIllegalArgumentException, MathIllegalStateException {
471 try {
472 TestProblem1 pb = new TestProblem1();
473 EmbeddedRungeKuttaIntegrator integrator = createIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
474 integrator.integrate(pb,
475 new ODEState(0.0, new double[pb.getDimension()]),
476 0.0);
477 fail("an exception should have been thrown");
478 } catch (MathIllegalArgumentException miae) {
479 assertEquals(LocalizedODEFormats.TOO_SMALL_INTEGRATION_INTERVAL,
480 miae.getSpecifier());
481 }
482 }
483
484 @Test
485 public abstract void testBackward();
486
487 protected void doTestBackward(final double epsilonLast, final double epsilonMaxValue,
488 final double epsilonMaxTime, final String name)
489 throws MathIllegalArgumentException, MathIllegalStateException {
490
491 TestProblem5 pb = new TestProblem5();
492 double minStep = 0;
493 double maxStep = FastMath.abs(pb.getFinalTime() - pb.getInitialState().getTime());
494 double scalAbsoluteTolerance = 1.0e-8;
495 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
496
497 EmbeddedRungeKuttaIntegrator integ = createIntegrator(minStep, maxStep,
498 scalAbsoluteTolerance,
499 scalRelativeTolerance);
500 TestProblemHandler handler = new TestProblemHandler(pb, integ);
501 integ.addStepHandler(handler);
502 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
503
504 assertEquals(0, handler.getLastError(), epsilonLast);
505 assertEquals(0, handler.getMaximalValueError(), epsilonMaxValue);
506 assertEquals(0, handler.getMaximalTimeError(), epsilonMaxTime);
507 assertEquals(name, integ.getName());
508
509 }
510
511 @Test
512 public abstract void testKepler();
513
514 protected void doTestKepler(double epsilon) {
515
516 final TestProblem3 pb = new TestProblem3(0.9);
517 double minStep = 1.0e-10;
518 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
519 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
520 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
521
522 AdaptiveStepsizeIntegrator integ = createIntegrator(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
523
524 try {
525 Method getStepSizeHelper = AdaptiveStepsizeIntegrator.class.getDeclaredMethod("getStepSizeHelper", (Class<?>[]) null);
526 getStepSizeHelper.setAccessible(true);
527 StepsizeHelper helper = (StepsizeHelper) getStepSizeHelper.invoke(integ, (Object[]) null);
528 integ.setInitialStepSize(-999);
529 assertEquals(-1.0, helper.getInitialStep(), 1.0e-10);
530 integ.setInitialStepSize(+999);
531 assertEquals(-1.0, helper.getInitialStep(), 1.0e-10);
532 } catch (NoSuchMethodException | IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
533 fail(e.getLocalizedMessage());
534 }
535
536 integ.addStepHandler(new KeplerHandler(pb, epsilon));
537 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
538 }
539
540 private static class KeplerHandler implements ODEStepHandler {
541 private double maxError;
542 private final TestProblem3 pb;
543 private final double epsilon;
544 public KeplerHandler(TestProblem3 pb, double epsilon) {
545 this.pb = pb;
546 this.epsilon = epsilon;
547 maxError = 0;
548 }
549 public void init(ODEStateAndDerivative state0, double t) {
550 maxError = 0;
551 }
552 public void handleStep(ODEStateInterpolator interpolator) {
553
554 ODEStateAndDerivative current = interpolator.getCurrentState();
555 double[] theoreticalY = pb.computeTheoreticalState(current.getTime());
556 double dx = current.getPrimaryState()[0] - theoreticalY[0];
557 double dy = current.getPrimaryState()[1] - theoreticalY[1];
558 double error = dx * dx + dy * dy;
559 if (error > maxError) {
560 maxError = error;
561 }
562 }
563 public void finish(ODEStateAndDerivative finalState) {
564 assertEquals(0.0, maxError, epsilon);
565 }
566 }
567
568 @Test
569 public abstract void testTorqueFreeMotionOmegaOnly();
570
571 protected void doTestTorqueFreeMotionOmegaOnly(double epsilon) {
572
573 final TestProblem7 pb = new TestProblem7();
574 double minStep = 1.0e-10;
575 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
576 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-8 };
577 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-10 };
578
579 EmbeddedRungeKuttaIntegrator integ = createIntegrator(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
580 integ.addStepHandler(new TorqueFreeOmegaOnlyHandler(pb, epsilon));
581 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
582 }
583
584 private static class TorqueFreeOmegaOnlyHandler implements ODEStepHandler {
585 private double maxError;
586 private final TestProblem7 pb;
587 private final double epsilon;
588 public TorqueFreeOmegaOnlyHandler(TestProblem7 pb, double epsilon) {
589 this.pb = pb;
590 this.epsilon = epsilon;
591 maxError = 0;
592 }
593 public void init(ODEStateAndDerivative state0, double t) {
594 maxError = 0;
595 }
596 public void handleStep(ODEStateInterpolator interpolator) {
597
598 ODEStateAndDerivative current = interpolator.getCurrentState();
599 double[] theoreticalY = pb.computeTheoreticalState(current.getTime());
600 double do1 = current.getPrimaryState()[0] - theoreticalY[0];
601 double do2 = current.getPrimaryState()[1] - theoreticalY[1];
602 double do3 = current.getPrimaryState()[2] - theoreticalY[2];
603 double error = do1 * do1 + do2 * do2 + do3 * do3;
604 if (error > maxError) {
605 maxError = error;
606 }
607 }
608 public void finish(ODEStateAndDerivative finalState) {
609 assertEquals(0.0, maxError, epsilon);
610 }
611 }
612
613 @Test
614
615
616
617 public abstract void testTorqueFreeMotion();
618
619 protected void doTestTorqueFreeMotion(double epsilonOmega, double epsilonQ) {
620
621 final double t0 = 0.0;
622 final double t1 = 20.0;
623 final Vector3D omegaBase = new Vector3D(5.0, 0.0, 4.0);
624 final Rotation rBase = new Rotation(0.9, 0.437, 0.0, 0.0, true);
625 final List<Double> inertiaBase = Arrays.asList(3.0 / 8.0, 1.0 / 2.0, 5.0 / 8.0);
626 double minStep = 1.0e-10;
627 double maxStep = t1 - t0;
628 double[] vecAbsoluteTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
629 double[] vecRelativeTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
630
631
632 Stream<TestProblem8> problems = Stream.empty();
633
634
635 problems = Stream.concat(problems,
636 permute(omegaBase).
637 map(omega -> new TestProblem8(t0, t1, omega, rBase,
638 inertiaBase.get(0), Vector3D.PLUS_I,
639 inertiaBase.get(1), Vector3D.PLUS_J,
640 inertiaBase.get(2), Vector3D.PLUS_K)));
641
642
643 problems = Stream.concat(problems,
644 permute(rBase).
645 map(r -> new TestProblem8(t0, t1, omegaBase, r,
646 inertiaBase.get(0), Vector3D.PLUS_I,
647 inertiaBase.get(1), Vector3D.PLUS_J,
648 inertiaBase.get(2), Vector3D.PLUS_K)));
649
650
651 problems = Stream.concat(problems,
652 permute(inertiaBase).
653 map(inertia -> new TestProblem8(t0, t1, omegaBase, rBase,
654 inertia.get(0), Vector3D.PLUS_I,
655 inertia.get(1), Vector3D.PLUS_J,
656 inertia.get(2), Vector3D.PLUS_K)));
657
658 problems.forEach(problem -> {
659 EmbeddedRungeKuttaIntegrator integ = createIntegrator(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
660 integ.addStepHandler(new TorqueFreeHandler(problem, epsilonOmega, epsilonQ));
661 integ.integrate(new ExpandableODE(problem), problem.getInitialState(), problem.getFinalTime());
662 });
663
664 }
665
666 @Test
667 public abstract void testTorqueFreeMotionIssue230();
668
669 protected void doTestTorqueFreeMotionIssue230(double epsilonOmega, double epsilonQ) {
670
671 double i1 = 3.0 / 8.0;
672 Vector3D a1 = Vector3D.PLUS_I;
673 double i2 = 5.0 / 8.0;
674 Vector3D a2 = Vector3D.PLUS_K;
675 double i3 = 1.0 / 2.0;
676 Vector3D a3 = Vector3D.PLUS_J;
677 Vector3D o0 = new Vector3D(5.0, 0.0, 4.0);
678 double o1 = Vector3D.dotProduct(o0, a1);
679 double o2 = Vector3D.dotProduct(o0, a2);
680 double o3 = Vector3D.dotProduct(o0, a3);
681 double e = 0.5 * (i1 * o1 * o1 + i2 * o2 * o2 + i3 * o3 * o3);
682 double r1 = FastMath.sqrt(2 * e * i1);
683 double r3 = FastMath.sqrt(2 * e * i3);
684 int n = 50;
685 for (int i = 0; i < n; ++i) {
686 SinCos sc = FastMath.sinCos(-0.5 * FastMath.PI * (i + 50) / 200);
687 Vector3D om = new Vector3D(r1 * sc.cos() / i1, a1, r3 * sc.sin() / i3, a3);
688 TestProblem8 problem = new TestProblem8(0, 20,
689 om,
690 new Rotation(0.9, 0.437, 0.0, 0.0, true),
691 i1, a1,
692 i2, a2,
693 i3, a3);
694
695 double minStep = 1.0e-10;
696 double maxStep = problem.getFinalTime() - problem.getInitialTime();
697 double[] vecAbsoluteTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
698 double[] vecRelativeTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
699
700 EmbeddedRungeKuttaIntegrator integ = createIntegrator(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
701 integ.addStepHandler(new TorqueFreeHandler(problem, epsilonOmega, epsilonQ));
702 integ.integrate(new ExpandableODE(problem), problem.getInitialState(), problem.getFinalTime() * 0.1);
703
704 }
705
706 }
707
708
709
710
711
712 private Stream<Vector3D> permute(final Vector3D v) {
713 return CombinatoricsUtils.
714 permutations(Arrays.asList(v.getX(), v.getY(), v.getZ())).
715 map(a -> new Vector3D(a.get(0), a.get(1), a.get(2)));
716 }
717
718
719
720
721
722 private Stream<Rotation> permute(final Rotation r) {
723 return CombinatoricsUtils.
724 permutations(Arrays.asList(r.getQ0(), r.getQ1(), r.getQ2(), r.getQ3())).
725 map(a -> new Rotation(a.get(0), a.get(1), a.get(2), a.get(3), false));
726 }
727
728
729
730
731
732 private Stream<List<Double>> permute(final List<Double> list) {
733 return CombinatoricsUtils.permutations(list);
734 }
735
736 private static class TorqueFreeHandler implements ODEStepHandler {
737 private double maxErrorOmega;
738 private double maxErrorQ;
739 private final TestProblem8 pb;
740 private final double epsilonOmega;
741 private final double epsilonQ;
742 private double outputStep;
743 private double current;
744
745 public TorqueFreeHandler(TestProblem8 pb, double epsilonOmega, double epsilonQ) {
746 this.pb = pb;
747 this.epsilonOmega = epsilonOmega;
748 this.epsilonQ = epsilonQ;
749 maxErrorOmega = 0;
750 maxErrorQ = 0;
751 outputStep = 0.01;
752 }
753
754 public void init(ODEStateAndDerivative state0, double t) {
755 maxErrorOmega = 0;
756 maxErrorQ = 0;
757 current = state0.getTime() - outputStep;
758 }
759
760 public void handleStep(ODEStateInterpolator interpolator) {
761
762 current += outputStep;
763 while (interpolator.getPreviousState().getTime() <= current &&
764 interpolator.getCurrentState().getTime() > current) {
765 ODEStateAndDerivative state = interpolator.getInterpolatedState(current);
766 final double[] theoretical = pb.computeTheoreticalState(state.getTime());
767 final double errorOmega = Vector3D.distance(new Vector3D(state.getPrimaryState()[0],
768 state.getPrimaryState()[1],
769 state.getPrimaryState()[2]),
770 new Vector3D(theoretical[0],
771 theoretical[1],
772 theoretical[2]));
773 maxErrorOmega = FastMath.max(maxErrorOmega, errorOmega);
774 final double errorQ = Rotation.distance(new Rotation(state.getPrimaryState()[3],
775 state.getPrimaryState()[4],
776 state.getPrimaryState()[5],
777 state.getPrimaryState()[6],
778 true),
779 new Rotation(theoretical[3],
780 theoretical[4],
781 theoretical[5],
782 theoretical[6],
783 true));
784 maxErrorQ = FastMath.max(maxErrorQ, errorQ);
785 current += outputStep;
786
787 }
788 }
789
790 public void finish(ODEStateAndDerivative finalState) {
791 assertEquals(0.0, maxErrorOmega, epsilonOmega);
792 assertEquals(0.0, maxErrorQ, epsilonQ);
793 }
794
795 }
796
797
798 @Test
799 public abstract void testMissedEndEvent();
800
801 protected void doTestMissedEndEvent(final double epsilonY, final double epsilonT) {
802 final double t0 = 1878250320.0000029;
803 final double tEvent = 1878250379.9999986;
804 final double[] k = { 1.0e-4, 1.0e-5, 1.0e-6 };
805 OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {
806
807 public int getDimension() {
808 return k.length;
809 }
810
811 public double[] computeDerivatives(double t, double[] y) {
812 final double[] yDot = new double[y.length];
813 for (int i = 0; i < y.length; ++i) {
814 yDot[i] = k[i] * y[i];
815 }
816 return yDot;
817 }
818 };
819
820 EmbeddedRungeKuttaIntegrator integrator = createIntegrator(0.0, 100.0, 1.0e-10, 1.0e-10);
821
822 double[] y0 = new double[k.length];
823 for (int i = 0; i < y0.length; ++i) {
824 y0[i] = i + 1;
825 }
826
827 integrator.setInitialStepSize(60.0);
828 ODEStateAndDerivative finalState = integrator.integrate(new ExpandableODE(ode), new ODEState(t0, y0), tEvent);
829 assertEquals(tEvent, finalState.getTime(), epsilonT);
830 double[] y = finalState.getPrimaryState();
831 for (int i = 0; i < y.length; ++i) {
832 assertEquals(y0[i] * FastMath.exp(k[i] * (finalState.getTime() - t0)), y[i], epsilonY);
833 }
834
835 integrator.setInitialStepSize(60.0);
836 integrator.addEventDetector(new ODEEventDetector() {
837 public AdaptableInterval getMaxCheckInterval() {
838 return (s, isForward)-> Double.POSITIVE_INFINITY;
839 }
840 public int getMaxIterationCount() {
841 return 100;
842 }
843 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
844 return new BracketingNthOrderBrentSolver(0, 1.0e-20, 0, 5);
845 }
846 public ODEEventHandler getHandler() {
847 return (state, detector, increasing) -> Action.CONTINUE;
848 }
849 public double g(ODEStateAndDerivative s) {
850 return s.getTime() - tEvent;
851 }
852 });
853 finalState = integrator.integrate(new ExpandableODE(ode), new ODEState(t0, y0), tEvent + 120);
854 assertEquals(tEvent + 120, finalState.getTime(), epsilonT);
855 y = finalState.getPrimaryState();
856 for (int i = 0; i < y.length; ++i) {
857 assertEquals(y0[i] * FastMath.exp(k[i] * (finalState.getTime() - t0)), y[i], epsilonY);
858 }
859
860 }
861
862 @Test
863 public void testTooLargeFirstStep() {
864
865 AdaptiveStepsizeIntegrator integ =
866 createIntegrator(0, Double.POSITIVE_INFINITY, Double.NaN, Double.NaN);
867 final double start = 0.0;
868 final double end = 0.001;
869 OrdinaryDifferentialEquation equations = new OrdinaryDifferentialEquation() {
870
871 public int getDimension() {
872 return 1;
873 }
874
875 public double[] computeDerivatives(double t, double[] y) {
876 assertTrue(t >= FastMath.nextAfter(start, Double.NEGATIVE_INFINITY));
877 assertTrue(t <= FastMath.nextAfter(end, Double.POSITIVE_INFINITY));
878 return new double[] { -100.0 * y[0] };
879 }
880
881 };
882
883 integ.setStepSizeControl(0, 1.0, 1.0e-6, 1.0e-8);
884 integ.integrate(equations, new ODEState(start, new double[] { 1.0 }), end);
885
886 }
887
888 @Test
889 public abstract void testVariableSteps();
890
891 protected void doTestVariableSteps(final double min, final double max) {
892
893 final TestProblem3 pb = new TestProblem3(0.9);
894 double minStep = 0;
895 double maxStep = pb.getFinalTime() - pb.getInitialTime();
896 double scalAbsoluteTolerance = 1.0e-8;
897 double scalRelativeTolerance = scalAbsoluteTolerance;
898
899 ODEIntegrator integ = createIntegrator(minStep, maxStep,
900 scalAbsoluteTolerance,
901 scalRelativeTolerance);
902 integ.addStepHandler(new VariableHandler(min, max));
903 double stopTime = integ.integrate(pb, pb.getInitialState(), pb.getFinalTime()).getTime();
904 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
905 }
906
907 @Test
908 public abstract void testUnstableDerivative();
909
910 protected void doTestUnstableDerivative(final double epsilon) {
911 final StepProblem stepProblem = new StepProblem((s, isForward) -> 999.0, 1.0e+12, 1000000, 0.0, 1.0, 2.0).
912 withMaxCheck(1.0).
913 withMaxIter(1000).
914 withThreshold(1.0e-12);
915 assertEquals(1.0, stepProblem.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
916 assertEquals(1000, stepProblem.getMaxIterationCount());
917 assertEquals(1.0e-12, stepProblem.getSolver().getAbsoluteAccuracy(), 1.0e-25);
918 assertNotNull(stepProblem.getHandler());
919 ODEIntegrator integ = createIntegrator(0.1, 10, 1.0e-12, 0.0);
920 integ.addEventDetector(stepProblem);
921 final ODEStateAndDerivative finalState =
922 integ.integrate(stepProblem, new ODEState(0.0, new double[] { 0.0 }), 10.0);
923 assertEquals(8.0, finalState.getPrimaryState()[0], epsilon);
924 }
925
926 @Test
927 public void testEventsScheduling() {
928
929 OrdinaryDifferentialEquation sincos = new OrdinaryDifferentialEquation() {
930
931 public int getDimension() {
932 return 2;
933 }
934
935 public double[] computeDerivatives(double t, double[] y) {
936 return new double[] { y[1], -y[0] };
937 }
938
939 };
940
941 SchedulingChecker sinChecker = new SchedulingChecker(0);
942 SchedulingChecker cosChecker = new SchedulingChecker(1);
943
944 ODEIntegrator integ = createIntegrator(0.001, 1.0, 1.0e-12, 0.0);
945 integ.addEventDetector(sinChecker);
946 integ.addStepHandler(sinChecker);
947 integ.addEventDetector(cosChecker);
948 integ.addStepHandler(cosChecker);
949 double t0 = 0.5;
950 double[] y0 = new double[] { FastMath.sin(t0), FastMath.cos(t0) };
951 double t = 10.0;
952 integ.integrate(sincos, new ODEState(t0, y0), t);
953
954 }
955
956 private static class SchedulingChecker implements ODEStepHandler, ODEEventDetector {
957
958 int index;
959 double tMin;
960
961 public SchedulingChecker(int index) {
962 this.index = index;
963 }
964
965 public AdaptableInterval getMaxCheckInterval() {
966 return (s, isForward)-> 0.01;
967 }
968
969 public int getMaxIterationCount() {
970 return 100;
971 }
972
973 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
974 return new BracketingNthOrderBrentSolver(0, 1.0e-7, 0, 5);
975 }
976
977 public void init(ODEStateAndDerivative s0, double t) {
978 tMin = s0.getTime();
979 }
980
981 public void handleStep(ODEStateInterpolator interpolator) {
982 tMin = interpolator.getCurrentState().getTime();
983 }
984
985 public double g(ODEStateAndDerivative s) {
986
987
988 assertTrue(s.getTime() >= tMin);
989 return s.getPrimaryState()[index];
990 }
991
992 public ODEEventHandler getHandler() {
993 return new ODEEventHandler() {
994 public Action eventOccurred(ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) {
995 return Action.RESET_STATE;
996 }
997
998 public ODEStateAndDerivative resetState(ODEEventDetector detector, ODEStateAndDerivative s) {
999
1000 return s;
1001 }
1002 };
1003 }
1004
1005 }
1006
1007 private static class VariableHandler implements ODEStepHandler {
1008 final double min;
1009 final double max;
1010 public VariableHandler(final double min, final double max) {
1011 this.min = min;
1012 this.max = max;
1013 firstTime = true;
1014 minStep = 0;
1015 maxStep = 0;
1016 }
1017 public void init(ODEStateAndDerivative s0, double t) {
1018 firstTime = true;
1019 minStep = 0;
1020 maxStep = 0;
1021 }
1022 public void handleStep(ODEStateInterpolator interpolator) {
1023
1024 double step = FastMath.abs(interpolator.getCurrentState().getTime() -
1025 interpolator.getPreviousState().getTime());
1026 if (firstTime) {
1027 minStep = FastMath.abs(step);
1028 maxStep = minStep;
1029 firstTime = false;
1030 } else {
1031 if (step < minStep) {
1032 minStep = step;
1033 }
1034 if (step > maxStep) {
1035 maxStep = step;
1036 }
1037 }
1038 }
1039 public void finish(ODEStateAndDerivative finalState) {
1040 assertEquals(min, minStep, 0.01 * min);
1041 assertEquals(max, maxStep, 0.01 * max);
1042 }
1043 private boolean firstTime = true;
1044 private double minStep = 0;
1045 private double maxStep = 0;
1046 }
1047
1048 @Test
1049 public void testWrongDerivative() {
1050 EmbeddedRungeKuttaIntegrator integrator =
1051 createIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
1052 OrdinaryDifferentialEquation equations =
1053 new OrdinaryDifferentialEquation() {
1054 public double[] computeDerivatives(double t, double[] y) {
1055 if (t < -0.5) {
1056 throw new LocalException();
1057 } else {
1058 throw new RuntimeException("oops");
1059 }
1060 }
1061 public int getDimension() {
1062 return 1;
1063 }
1064 };
1065
1066 try {
1067 integrator.integrate(equations, new ODEState(-1.0, new double[1]), 0.0);
1068 fail("an exception should have been thrown");
1069 } catch(LocalException de) {
1070
1071 }
1072
1073 try {
1074 integrator.integrate(equations, new ODEState(0.0, new double[1]), 1.0);
1075 fail("an exception should have been thrown");
1076 } catch(RuntimeException de) {
1077
1078 }
1079
1080 }
1081
1082 @Test
1083 public abstract void testPartialDerivatives();
1084
1085 protected void doTestPartialDerivatives(final double epsilonY,
1086 final double epsilonPartials) {
1087
1088 double omega = 1.3;
1089 double t0 = 1.3;
1090 double[] y0 = new double[] { 3.0, 4.0 };
1091 double t = 6.0;
1092 SinCosJacobiansProvider sinCos = new SinCosJacobiansProvider(omega);
1093
1094 ExpandableODE expandable = new ExpandableODE(sinCos);
1095 VariationalEquation ve = new VariationalEquation(expandable, sinCos);
1096 ODEState initialState = ve.setUpInitialState(new ODEState(t0, y0));
1097
1098 EmbeddedRungeKuttaIntegrator integrator =
1099 createIntegrator(0.001 * (t - t0), t - t0, 1.0e-12, 1.0e-12);
1100 ODEStateAndDerivative finalState = integrator.integrate(expandable, initialState, t);
1101
1102
1103 int n = sinCos.getDimension();
1104 int p = sinCos.getParametersNames().size();
1105 assertEquals(n, finalState.getPrimaryStateDimension());
1106 assertEquals(n + n * (n + p), finalState.getCompleteStateDimension());
1107
1108
1109 for (int i = 0; i < sinCos.getDimension(); ++i) {
1110 assertEquals(sinCos.theoreticalY(t)[i], finalState.getPrimaryState()[i], epsilonY);
1111 }
1112
1113
1114 double[][] dydy0 = ve.extractMainSetJacobian(finalState);
1115 for (int i = 0; i < dydy0.length; ++i) {
1116 for (int j = 0; j < dydy0[i].length; ++j) {
1117 assertEquals(sinCos.exactDyDy0(t)[i][j], dydy0[i][j], epsilonPartials);
1118 }
1119 }
1120
1121 double[] dydp0 = ve.extractParameterJacobian(finalState, SinCosJacobiansProvider.OMEGA_PARAMETER);
1122 for (int i = 0; i < dydp0.length; ++i) {
1123 assertEquals(sinCos.exactDyDomega(t)[i], dydp0[i], epsilonPartials);
1124 }
1125
1126 }
1127
1128 @Test
1129 public abstract void testSecondaryEquations();
1130
1131 protected void doTestSecondaryEquations(final double epsilonSinCos,
1132 final double epsilonLinear) {
1133 OrdinaryDifferentialEquation sinCos = new OrdinaryDifferentialEquation() {
1134
1135 @Override
1136 public int getDimension() {
1137 return 2;
1138 }
1139
1140 @Override
1141 public double[] computeDerivatives(double t, double[] y) {
1142 return new double[] { y[1], -y[0] };
1143 }
1144
1145 };
1146
1147 SecondaryODE linear = new SecondaryODE() {
1148
1149 @Override
1150 public int getDimension() {
1151 return 1;
1152 }
1153
1154 @Override
1155 public double[] computeDerivatives(double t, double[] primary, double[] primaryDot, double[] secondary) {
1156 return new double[] { -1 };
1157 }
1158
1159 };
1160
1161 ExpandableODE expandable = new ExpandableODE(sinCos);
1162 expandable.addSecondaryEquations(linear);
1163
1164 ODEIntegrator integrator = createIntegrator(0.001, 1.0, 1.0e-12, 1.0e-12);
1165 final double[] max = new double[2];
1166 integrator.addStepHandler(new ODEStepHandler() {
1167 @Override
1168 public void handleStep(ODEStateInterpolator interpolator) {
1169 for (int i = 0; i <= 10; ++i) {
1170 double tPrev = interpolator.getPreviousState().getTime();
1171 double tCurr = interpolator.getCurrentState().getTime();
1172 double t = (tPrev * (10 - i) + tCurr * i) / 10;
1173 ODEStateAndDerivative state = interpolator.getInterpolatedState(t);
1174 assertEquals(2, state.getPrimaryStateDimension());
1175 assertEquals(1, state.getNumberOfSecondaryStates());
1176 assertEquals(2, state.getSecondaryStateDimension(0));
1177 assertEquals(1, state.getSecondaryStateDimension(1));
1178 assertEquals(3, state.getCompleteStateDimension());
1179 max[0] = FastMath.max(max[0],
1180 FastMath.abs(FastMath.sin(t) - state.getPrimaryState()[0]));
1181 max[0] = FastMath.max(max[0],
1182 FastMath.abs(FastMath.cos(t) - state.getPrimaryState()[1]));
1183 max[1] = FastMath.max(max[1],
1184 FastMath.abs(1 - t - state.getSecondaryState(1)[0]));
1185 }
1186 }
1187 });
1188
1189 double[] primary0 = new double[] { 0.0, 1.0 };
1190 double[][] secondary0 = new double[][] { { 1.0 } };
1191 ODEState initialState = new ODEState(0.0, primary0, secondary0);
1192
1193 ODEStateAndDerivative finalState =
1194 integrator.integrate(expandable, initialState, 10.0);
1195 assertEquals(10.0, finalState.getTime(), 1.0e-12);
1196 assertEquals(0, max[0], epsilonSinCos);
1197 assertEquals(0, max[1], epsilonLinear);
1198
1199 }
1200
1201 @Test
1202 public void testNaNAppearing() {
1203 try {
1204 ODEIntegrator integ = createIntegrator(0.01, 1.0, 0.1, 0.1);
1205 integ.integrate(new OrdinaryDifferentialEquation() {
1206 public int getDimension() {
1207 return 1;
1208 }
1209 public double[] computeDerivatives(double t, double[] y) {
1210 return new double[] { FastMath.log(t) };
1211 }
1212 }, new ODEState(1.0, new double[] { 1.0 }), -1.0);
1213 fail("an exception should have been thrown");
1214 } catch (MathIllegalStateException mise) {
1215 assertEquals(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION, mise.getSpecifier());
1216 assertTrue(((Double) mise.getParts()[0]).doubleValue() <= 0.0);
1217 }
1218 }
1219
1220 @Test
1221 public void testInfiniteIntegration() {
1222 ODEIntegrator integ = createIntegrator(0.01, 1.0, 0.1, 0.1);
1223 TestProblem1 pb = new TestProblem1();
1224 double convergence = 1e-6;
1225 integ.addEventDetector(new ODEEventDetector() {
1226 @Override
1227 public AdaptableInterval getMaxCheckInterval() {
1228 return (s, isForward)-> Double.POSITIVE_INFINITY;
1229 }
1230 @Override
1231 public int getMaxIterationCount() {
1232 return 1000;
1233 }
1234 @Override
1235 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
1236 return new BracketingNthOrderBrentSolver(0, convergence, 0, 5);
1237 }
1238 @Override
1239 public double g(ODEStateAndDerivative state) {
1240 return state.getTime() - pb.getFinalTime();
1241 }
1242 @Override
1243 public ODEEventHandler getHandler() {
1244 return (state, detector, increasing) -> Action.STOP;
1245 }
1246 });
1247 ODEStateAndDerivative finalState = integ.integrate(pb, pb.getInitialState(), Double.POSITIVE_INFINITY);
1248 assertEquals(pb.getFinalTime(), finalState.getTime(), convergence);
1249 }
1250
1251 private static class SinCosJacobiansProvider implements ODEJacobiansProvider {
1252
1253 public static String OMEGA_PARAMETER = "omega";
1254
1255 private final double omega;
1256 private double r;
1257 private double alpha;
1258
1259 private double dRdY00;
1260 private double dRdY01;
1261 private double dAlphadOmega;
1262 private double dAlphadY00;
1263 private double dAlphadY01;
1264
1265 protected SinCosJacobiansProvider(final double omega) {
1266 this.omega = omega;
1267 }
1268
1269 public int getDimension() {
1270 return 2;
1271 }
1272
1273 public void init(final double t0, final double[] y0,
1274 final double finalTime) {
1275
1276
1277
1278 final double r2 = y0[0] * y0[0] + y0[1] * y0[1];
1279
1280 this.r = FastMath.sqrt(r2);
1281 this.dRdY00 = y0[0] / r;
1282 this.dRdY01 = y0[1] / r;
1283
1284 this.alpha = FastMath.atan2(y0[0], y0[1]) - t0 * omega;
1285 this.dAlphadOmega = -t0;
1286 this.dAlphadY00 = y0[1] / r2;
1287 this.dAlphadY01 = -y0[0] / r2;
1288
1289 }
1290
1291 @Override
1292 public double[] computeDerivatives(final double t, final double[] y) {
1293 return new double[] {
1294 omega * y[1],
1295 omega * -y[0]
1296 };
1297 }
1298
1299 @Override
1300 public double[][] computeMainStateJacobian(final double t, final double[] y, final double[] yDot) {
1301
1302 return new double[][] {
1303 { 0.0, omega },
1304 { -omega, 0.0 }
1305 };
1306 }
1307
1308 @Override
1309 public double[] computeParameterJacobian(final double t, final double[] y, final double[] yDot,
1310 final String paramName)
1311 throws MathIllegalArgumentException, MathIllegalStateException {
1312 if (!isSupported(paramName)) {
1313 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, paramName);
1314 }
1315
1316 return new double[] {
1317 y[1],
1318 -y[0]
1319 };
1320 }
1321
1322 public double[] theoreticalY(final double t) {
1323 final double theta = omega * t + alpha;
1324 return new double[] {
1325 r * FastMath.sin(theta), r * FastMath.cos(theta)
1326 };
1327 }
1328
1329 public double[][] exactDyDy0(final double t) {
1330
1331
1332 final double theta = omega * t + alpha;
1333 final double sin = FastMath.sin(theta);
1334 final double cos = FastMath.cos(theta);
1335 final double y0 = r * sin;
1336 final double y1 = r * cos;
1337
1338 return new double[][] {
1339 { dRdY00 * sin + y1 * dAlphadY00, dRdY01 * sin + y1 * dAlphadY01 },
1340 { dRdY00 * cos - y0 * dAlphadY00, dRdY01 * cos - y0 * dAlphadY01 }
1341 };
1342
1343 }
1344
1345 public double[] exactDyDomega(final double t) {
1346
1347
1348 final double theta = omega * t + alpha;
1349 final double sin = FastMath.sin(theta);
1350 final double cos = FastMath.cos(theta);
1351 final double y0 = r * sin;
1352 final double y1 = r * cos;
1353
1354 return new double[] {
1355 y1 * (t + dAlphadOmega),
1356 -y0 * (t + dAlphadOmega)
1357 };
1358
1359 }
1360
1361 @Override
1362 public List<String> getParametersNames() {
1363 return Arrays.asList(OMEGA_PARAMETER);
1364 }
1365
1366 @Override
1367 public boolean isSupported(final String name) {
1368 return OMEGA_PARAMETER.equals(name);
1369 }
1370
1371 }
1372
1373 }