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