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.Constructor;
22 import java.lang.reflect.InvocationTargetException;
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.CalculusFieldElement;
29 import org.hipparchus.Field;
30 import org.hipparchus.analysis.differentiation.DSFactory;
31 import org.hipparchus.analysis.differentiation.DerivativeStructure;
32 import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
33 import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
34 import org.hipparchus.exception.MathIllegalArgumentException;
35 import org.hipparchus.exception.MathIllegalStateException;
36 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
37 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
38 import org.hipparchus.ode.FieldExpandableODE;
39 import org.hipparchus.ode.FieldODEIntegrator;
40 import org.hipparchus.ode.FieldODEState;
41 import org.hipparchus.ode.FieldODEStateAndDerivative;
42 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
43 import org.hipparchus.ode.FieldSecondaryODE;
44 import org.hipparchus.ode.TestFieldProblem1;
45 import org.hipparchus.ode.TestFieldProblem3;
46 import org.hipparchus.ode.TestFieldProblem4;
47 import org.hipparchus.ode.TestFieldProblem5;
48 import org.hipparchus.ode.TestFieldProblem7;
49 import org.hipparchus.ode.TestFieldProblem8;
50 import org.hipparchus.ode.TestFieldProblemHandler;
51 import org.hipparchus.ode.events.Action;
52 import org.hipparchus.ode.events.FieldAdaptableInterval;
53 import org.hipparchus.ode.events.FieldODEEventDetector;
54 import org.hipparchus.ode.events.FieldODEEventHandler;
55 import org.hipparchus.ode.events.FieldODEStepEndHandler;
56 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
57 import org.hipparchus.ode.sampling.FieldODEStepHandler;
58 import org.hipparchus.util.Binary64;
59 import org.hipparchus.util.Binary64Field;
60 import org.hipparchus.util.CombinatoricsUtils;
61 import org.hipparchus.util.FastMath;
62 import org.hipparchus.util.MathArrays;
63 import org.hipparchus.util.SinCos;
64 import org.junit.Assert;
65 import org.junit.Test;
66
67 public abstract class EmbeddedRungeKuttaFieldIntegratorAbstractTest {
68
69 protected abstract <T extends CalculusFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
70 createIntegrator(Field<T> field, final double minStep, final double maxStep,
71 final double scalAbsoluteTolerance, final double scalRelativeTolerance);
72
73 protected abstract <T extends CalculusFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
74 createIntegrator(Field<T> field, final double minStep, final double maxStep,
75 final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
76
77 @Test
78 public abstract void testNonFieldIntegratorConsistency();
79
80 protected <T extends CalculusFieldElement<T>> void doTestNonFieldIntegratorConsistency(final Field<T> field) {
81 try {
82
83
84 EmbeddedRungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, 0.001, 1.0, 1.0, 1.0);
85 T[][] fieldA = fieldIntegrator.getA();
86 T[] fieldB = fieldIntegrator.getB();
87 T[] fieldC = fieldIntegrator.getC();
88
89 String fieldName = fieldIntegrator.getClass().getName();
90 String regularName = fieldName.replaceAll("Field", "");
91
92
93 @SuppressWarnings("unchecked")
94 Constructor<EmbeddedRungeKuttaIntegrator> constructor =
95 (Constructor<EmbeddedRungeKuttaIntegrator>) Class.forName(regularName).getConstructor(Double.TYPE,
96 Double.TYPE,
97 Double.TYPE,
98 Double.TYPE);
99 final EmbeddedRungeKuttaIntegrator regularIntegrator =
100 constructor.newInstance(0.0, 1.0, 1.0e-10, 1.0e-10);
101 double[][] regularA = regularIntegrator.getA();
102 double[] regularB = regularIntegrator.getB();
103 double[] regularC = regularIntegrator.getC();
104
105 Assert.assertEquals(regularA.length, fieldA.length);
106 for (int i = 0; i < regularA.length; ++i) {
107 checkArray(regularA[i], fieldA[i]);
108 }
109 checkArray(regularB, fieldB);
110 checkArray(regularC, fieldC);
111
112 } catch (ClassNotFoundException | IllegalAccessException | IllegalArgumentException |
113 SecurityException | NoSuchMethodException | InvocationTargetException |
114 InstantiationException e) {
115 Assert.fail(e.getLocalizedMessage());
116 }
117 }
118
119 private <T extends CalculusFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) {
120 Assert.assertEquals(regularArray.length, fieldArray.length);
121 for (int i = 0; i < regularArray.length; ++i) {
122 if (regularArray[i] == 0) {
123 Assert.assertTrue(0.0 == fieldArray[i].getReal());
124 } else {
125 Assert.assertEquals(regularArray[i], fieldArray[i].getReal(), FastMath.ulp(regularArray[i]));
126 }
127 }
128 }
129
130 @Test
131 public abstract void testForwardBackwardExceptions();
132
133 protected <T extends CalculusFieldElement<T>> void doTestForwardBackwardExceptions(final Field<T> field) {
134 FieldOrdinaryDifferentialEquation<T> equations = new FieldOrdinaryDifferentialEquation<T>() {
135
136 public int getDimension() {
137 return 1;
138 }
139
140 public void init(T t0, T[] y0, T t) {
141 }
142
143 public T[] computeDerivatives(T t, T[] y) {
144 if (t.getReal() < -0.5) {
145 throw new LocalException();
146 } else {
147 throw new RuntimeException("oops");
148 }
149 }
150 };
151
152 EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0.0, 1.0, 1.0e-10, 1.0e-10);
153
154 try {
155 integrator.integrate(new FieldExpandableODE<T>(equations),
156 new FieldODEState<T>(field.getOne().negate(),
157 MathArrays.buildArray(field, 1)),
158 field.getZero());
159 Assert.fail("an exception should have been thrown");
160 } catch(LocalException de) {
161
162 }
163
164 try {
165 integrator.integrate(new FieldExpandableODE<T>(equations),
166 new FieldODEState<T>(field.getZero(),
167 MathArrays.buildArray(field, 1)),
168 field.getOne());
169 Assert.fail("an exception should have been thrown");
170 } catch(RuntimeException de) {
171
172 }
173 }
174
175 protected static class LocalException extends RuntimeException {
176 private static final long serialVersionUID = 20151208L;
177 }
178
179 @Test(expected=MathIllegalArgumentException.class)
180 public abstract void testMinStep();
181
182 protected <T extends CalculusFieldElement<T>> void doTestMinStep(final Field<T> field)
183 throws MathIllegalArgumentException {
184
185 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
186 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
187 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
188 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
189 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
190
191 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
192 vecAbsoluteTolerance, vecRelativeTolerance);
193 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
194 integ.addStepHandler(handler);
195 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
196 Assert.fail("an exception should have been thrown");
197
198 }
199
200 @Test
201 public abstract void testIncreasingTolerance();
202
203 protected <T extends CalculusFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
204 double factor,
205 double epsilon) {
206
207 int previousCalls = Integer.MAX_VALUE;
208 for (int i = -12; i < -2; ++i) {
209 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
210 double minStep = 0;
211 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
212 double scalAbsoluteTolerance = FastMath.pow(10.0, i);
213 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
214
215 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
216 scalAbsoluteTolerance, scalRelativeTolerance);
217 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
218 integ.addStepHandler(handler);
219 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
220
221 Assert.assertTrue(handler.getMaximalValueError().getReal() < (factor * scalAbsoluteTolerance));
222 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilon);
223
224 int calls = pb.getCalls();
225 Assert.assertEquals(integ.getEvaluations(), calls);
226 Assert.assertTrue(calls <= previousCalls);
227 previousCalls = calls;
228
229 }
230
231 }
232
233 @Test
234 public abstract void testEvents();
235
236 protected <T extends CalculusFieldElement<T>> void doTestEvents(final Field<T> field,
237 final double epsilonMaxValue,
238 final String name) {
239
240 TestFieldProblem4<T> pb = new TestFieldProblem4<T>(field);
241 double minStep = 0;
242 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
243 double scalAbsoluteTolerance = 1.0e-8;
244 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
245
246 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
247 scalAbsoluteTolerance, scalRelativeTolerance);
248 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
249 integ.addStepHandler(handler);
250 double convergence = 1.0e-8 * maxStep;
251 FieldODEEventDetector<T>[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY,
252 field.getZero().newInstance(convergence),
253 1000);
254 for (int l = 0; l < functions.length; ++l) {
255 integ.addEventDetector(functions[l]);
256 }
257 List<FieldODEEventDetector<T>> detectors = new ArrayList<>(integ.getEventDetectors());
258 Assert.assertEquals(functions.length, integ.getEventDetectors().size());
259
260 for (int i = 0; i < detectors.size(); ++i) {
261 Assert.assertSame(functions[i], detectors.get(i).getHandler());
262 Assert.assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null), 1.0);
263 Assert.assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy().getReal(), 1.0e-15 * convergence);
264 Assert.assertEquals(1000, detectors.get(i).getMaxIterationCount());
265 }
266
267 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
268
269 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
270 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), convergence);
271 Assert.assertEquals(12.0, handler.getLastTime().getReal(), convergence);
272 Assert.assertEquals(name, integ.getName());
273 integ.clearEventDetectors();
274 Assert.assertEquals(0, integ.getEventDetectors().size());
275
276 }
277
278 @Test
279 public abstract void testStepEnd();
280
281 protected <T extends CalculusFieldElement<T>> void doTestStepEnd(final Field<T> field,
282 final int expectedCount,
283 final String name) {
284
285 TestFieldProblem4<T> pb = new TestFieldProblem4<T>(field);
286 double minStep = 0;
287 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
288 double scalAbsoluteTolerance = 1.0e-8;
289 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
290
291 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
292 scalAbsoluteTolerance, scalRelativeTolerance);
293 double convergence = 1.0e-8 * maxStep;
294 FieldODEEventDetector<T>[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY,
295 field.getZero().newInstance(convergence),
296 1000);
297 for (int l = 0; l < functions.length; ++l) {
298 integ.addEventDetector(functions[l]);
299 }
300 List<FieldODEEventDetector<T>> detectors = new ArrayList<>(integ.getEventDetectors());
301 Assert.assertEquals(functions.length, integ.getEventDetectors().size());
302
303 for (int i = 0; i < detectors.size(); ++i) {
304 Assert.assertSame(functions[i], detectors.get(i).getHandler());
305 Assert.assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null), 1.0);
306 Assert.assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy().getReal(), 1.0e-15 * convergence);
307 Assert.assertEquals(1000, detectors.get(i).getMaxIterationCount());
308 }
309
310 final StepCounter<T> counter = new StepCounter<>(expectedCount + 10, Action.STOP);
311 integ.addStepEndHandler(counter);
312 Assert.assertEquals(1, integ.getStepEndHandlers().size());
313 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
314
315 Assert.assertEquals(expectedCount, counter.count);
316 Assert.assertEquals(name, integ.getName());
317 integ.clearEventDetectors();
318 Assert.assertEquals(0, integ.getEventDetectors().size());
319 integ.clearStepEndHandlers();
320 Assert.assertEquals(0, integ.getStepEndHandlers().size());
321
322 }
323
324 @Test
325 public abstract void testStopAfterStep();
326
327 protected <T extends CalculusFieldElement<T>> void doTestStopAfterStep(final Field<T> field,
328 final int count,
329 final double expectedTime) {
330
331 TestFieldProblem4<T> pb = new TestFieldProblem4<T>(field);
332 double minStep = 0;
333 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
334 double scalAbsoluteTolerance = 1.0e-8;
335 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
336
337 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
338 scalAbsoluteTolerance, scalRelativeTolerance);
339 double convergence = 1.0e-8 * maxStep;
340 FieldODEEventDetector<T>[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY,
341 field.getZero().newInstance(convergence),
342 1000);
343 for (int l = 0; l < functions.length; ++l) {
344 integ.addEventDetector(functions[l]);
345 }
346 List<FieldODEEventDetector<T>> detectors = new ArrayList<>(integ.getEventDetectors());
347 Assert.assertEquals(functions.length, integ.getEventDetectors().size());
348
349 for (int i = 0; i < detectors.size(); ++i) {
350 Assert.assertSame(functions[i], detectors.get(i).getHandler());
351 Assert.assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null), 1.0);
352 Assert.assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy().getReal(), 1.0e-15 * convergence);
353 Assert.assertEquals(1000, detectors.get(i).getMaxIterationCount());
354 }
355
356 final StepCounter<T> counter = new StepCounter<>(count, Action.STOP);
357 integ.addStepEndHandler(counter);
358 Assert.assertEquals(1, integ.getStepEndHandlers().size());
359 FieldODEStateAndDerivative<T> finalState = integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
360
361 Assert.assertEquals(count, counter.count);
362 Assert.assertEquals(expectedTime, finalState.getTime().getReal(), 1.0e-6);
363
364 }
365
366 @Test
367 public abstract void testResetAfterStep();
368
369 protected <T extends CalculusFieldElement<T>> void doTestResetAfterStep(final Field<T> field,
370 final int resetCount,
371 final int expectedCount) {
372
373 TestFieldProblem4<T> pb = new TestFieldProblem4<T>(field);
374 double minStep = 0;
375 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
376 double scalAbsoluteTolerance = 1.0e-8;
377 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
378
379 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
380 scalAbsoluteTolerance, scalRelativeTolerance);
381 double convergence = 1.0e-8 * maxStep;
382 FieldODEEventDetector<T>[] functions = pb.getEventDetectors(Double.POSITIVE_INFINITY,
383 field.getZero().newInstance(convergence),
384 1000);
385 for (int l = 0; l < functions.length; ++l) {
386 integ.addEventDetector(functions[l]);
387 }
388 List<FieldODEEventDetector<T>> detectors = new ArrayList<>(integ.getEventDetectors());
389 Assert.assertEquals(functions.length, integ.getEventDetectors().size());
390
391 for (int i = 0; i < detectors.size(); ++i) {
392 Assert.assertSame(functions[i], detectors.get(i).getHandler());
393 Assert.assertEquals(Double.POSITIVE_INFINITY, detectors.get(i).getMaxCheckInterval().currentInterval(null), 1.0);
394 Assert.assertEquals(convergence, detectors.get(i).getSolver().getAbsoluteAccuracy().getReal(), 1.0e-15 * convergence);
395 Assert.assertEquals(1000, detectors.get(i).getMaxIterationCount());
396 }
397
398 final StepCounter<T> counter = new StepCounter<>(resetCount, Action.RESET_STATE);
399 integ.addStepEndHandler(counter);
400 Assert.assertEquals(1, integ.getStepEndHandlers().size());
401 FieldODEStateAndDerivative<T> finalState = integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
402
403 Assert.assertEquals(expectedCount, counter.count);
404 Assert.assertEquals(12.0, finalState.getTime().getReal(), 1.0e-6);
405 for (int i = 0; i < finalState.getPrimaryStateDimension(); ++i) {
406 Assert.assertEquals(0.0, finalState.getPrimaryState()[i].getReal(), 1.0e-15);
407 Assert.assertEquals(0.0, finalState.getPrimaryDerivative()[i].getReal(), 1.0e-15);
408 }
409
410 }
411
412 private static class StepCounter<T extends CalculusFieldElement<T>> implements FieldODEStepEndHandler<T> {
413 final int max;
414 final Action actionAtMax;
415 int count;
416 StepCounter(final int max, final Action actionAtMax) {
417 this.max = max;
418 this.actionAtMax = actionAtMax;
419 this.count = 0;
420 }
421 public Action stepEndOccurred(final FieldODEStateAndDerivative<T> state, final boolean forward) {
422 return ++count == max ? actionAtMax : Action.CONTINUE;
423 }
424 public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
425 return new FieldODEState<>(state.getTime(),
426 MathArrays.buildArray(state.getTime().getField(),
427 state.getPrimaryStateDimension()));
428 }
429 }
430
431 @Test(expected=LocalException.class)
432 public abstract void testEventsErrors();
433
434 protected <T extends CalculusFieldElement<T>> void doTestEventsErrors(final Field<T> field)
435 throws LocalException {
436 final TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
437 double minStep = 0;
438 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
439 double scalAbsoluteTolerance = 1.0e-8;
440 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
441
442 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
443 scalAbsoluteTolerance, scalRelativeTolerance);
444 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
445 integ.addStepHandler(handler);
446
447 integ.addEventDetector(new FieldODEEventDetector<T>() {
448 public FieldAdaptableInterval<T> getMaxCheckInterval() {
449 return s -> Double.POSITIVE_INFINITY;
450 }
451 public int getMaxIterationCount() {
452 return 1000;
453 }
454 public BracketedRealFieldUnivariateSolver<T> getSolver() {
455 return new FieldBracketingNthOrderBrentSolver<T>(field.getZero(),
456 field.getZero().newInstance(1.0e-8 * maxStep),
457 field.getZero(),
458 5);
459 }
460 public FieldODEEventHandler<T> getHandler() {
461 return (state, detector, increasing) -> Action.CONTINUE;
462 }
463 public T g(FieldODEStateAndDerivative<T> state) {
464 T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
465 T offset = state.getTime().subtract(middle);
466 if (offset.getReal() > 0) {
467 throw new LocalException();
468 }
469 return offset;
470 }
471 });
472
473 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
474
475 }
476
477 @Test
478 public abstract void testEventsNoConvergence();
479
480 protected <T extends CalculusFieldElement<T>> void doTestEventsNoConvergence(final Field<T> field){
481
482 final TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
483 double minStep = 0;
484 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
485 double scalAbsoluteTolerance = 1.0e-8;
486 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
487
488 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
489 scalAbsoluteTolerance, scalRelativeTolerance);
490 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
491 integ.addStepHandler(handler);
492
493 integ.addEventDetector(new FieldODEEventDetector<T>() {
494 public FieldAdaptableInterval<T> getMaxCheckInterval() {
495 return s -> Double.POSITIVE_INFINITY;
496 }
497 public int getMaxIterationCount() {
498 return 3;
499 }
500 public BracketedRealFieldUnivariateSolver<T> getSolver() {
501 return new FieldBracketingNthOrderBrentSolver<T>(field.getZero(),
502 field.getZero().newInstance(1.0e-8 * maxStep),
503 field.getZero(),
504 5);
505 }
506 public FieldODEEventHandler<T> getHandler() {
507 return (state, detector, increasing) -> Action.CONTINUE;
508 }
509 public T g(FieldODEStateAndDerivative<T> state) {
510 T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
511 T offset = state.getTime().subtract(middle);
512 return (offset.getReal() > 0) ? offset.add(0.5) : offset.subtract(0.5);
513 }
514 });
515
516 try {
517 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
518 Assert.fail("an exception should have been thrown");
519 } catch (MathIllegalStateException mcee) {
520
521 }
522
523 }
524
525 @Test
526 public abstract void testSanityChecks();
527
528 protected <T extends CalculusFieldElement<T>> void doTestSanityChecks(Field<T> field) {
529 TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field);
530 try {
531 EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0,
532 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
533 new double[4], new double[4]);
534 integrator.integrate(new FieldExpandableODE<T>(pb),
535 new FieldODEState<T>(pb.getInitialState().getTime(),
536 MathArrays.buildArray(field, 6)),
537 pb.getFinalTime());
538 Assert.fail("an exception should have been thrown");
539 } catch(MathIllegalArgumentException ie) {
540 }
541 try {
542 EmbeddedRungeKuttaFieldIntegrator<T> integrator =
543 createIntegrator(field, 0,
544 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
545 new double[2], new double[4]);
546 integrator.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
547 Assert.fail("an exception should have been thrown");
548 } catch(MathIllegalArgumentException ie) {
549 }
550 try {
551 EmbeddedRungeKuttaFieldIntegrator<T> integrator =
552 createIntegrator(field, 0,
553 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
554 new double[4], new double[4]);
555 integrator.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getInitialState().getTime());
556 Assert.fail("an exception should have been thrown");
557 } catch(MathIllegalArgumentException ie) {
558 }
559 }
560
561 @Test
562 public abstract void testBackward();
563
564 protected <T extends CalculusFieldElement<T>> void doTestBackward(Field<T> field,
565 final double epsilonLast,
566 final double epsilonMaxValue,
567 final double epsilonMaxTime,
568 final String name)
569 throws MathIllegalArgumentException, MathIllegalStateException {
570
571 TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
572 double minStep = 0;
573 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).norm();
574 double scalAbsoluteTolerance = 1.0e-8;
575 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
576
577 EmbeddedRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
578 scalAbsoluteTolerance,
579 scalRelativeTolerance);
580 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
581 integ.addStepHandler(handler);
582 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
583
584 Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast);
585 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
586 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
587 Assert.assertEquals(name, integ.getName());
588
589 }
590
591 @Test
592 public abstract void testKepler();
593
594 protected <T extends CalculusFieldElement<T>> void doTestKepler(Field<T> field, double epsilon) {
595
596 final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field.getZero().add(0.9));
597 double minStep = 0;
598 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
599 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
600 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
601
602 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
603 vecAbsoluteTolerance, vecRelativeTolerance);
604 integ.addStepHandler(new KeplerHandler<T>(pb, epsilon));
605 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
606 }
607
608 private static class KeplerHandler<T extends CalculusFieldElement<T>> implements FieldODEStepHandler<T> {
609 private T maxError;
610 private final TestFieldProblem3<T> pb;
611 private final double epsilon;
612 public KeplerHandler(TestFieldProblem3<T> pb, double epsilon) {
613 this.pb = pb;
614 this.epsilon = epsilon;
615 maxError = pb.getField().getZero();
616 }
617 public void init(FieldODEStateAndDerivative<T> state0, T t) {
618 maxError = pb.getField().getZero();
619 }
620 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
621
622 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
623 T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
624 T dx = current.getPrimaryState()[0].subtract(theoreticalY[0]);
625 T dy = current.getPrimaryState()[1].subtract(theoreticalY[1]);
626 T error = dx.multiply(dx).add(dy.multiply(dy));
627 if (error.subtract(maxError).getReal() > 0) {
628 maxError = error;
629 }
630 }
631 public void finish(FieldODEStateAndDerivative<T> finalState) {
632 Assert.assertEquals(0.0, maxError.getReal(), epsilon);
633 }
634 }
635
636 @Test
637 public abstract void testTorqueFreeMotionOmegaOnly();
638
639 protected <T extends CalculusFieldElement<T>> void doTestTorqueFreeMotionOmegaOnly(Field<T> field, double epsilon) {
640
641 final TestFieldProblem7<T> pb = new TestFieldProblem7<>(field);
642 double minStep = 1.0e-10;
643 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
644 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-8 };
645 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-10 };
646
647 FieldODEIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
648 vecAbsoluteTolerance, vecRelativeTolerance);
649 integ.addStepHandler(new TorqueFreeHandlerOmegaOnly<>(pb, epsilon));
650 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
651 }
652
653 private static class TorqueFreeHandlerOmegaOnly<T extends CalculusFieldElement<T>> implements FieldODEStepHandler<T> {
654 private T maxError;
655 private final TestFieldProblem7<T> pb;
656 private final double epsilon;
657 public TorqueFreeHandlerOmegaOnly(TestFieldProblem7<T> pb, double epsilon) {
658 this.pb = pb;
659 this.epsilon = epsilon;
660 maxError = pb.getField().getZero();
661 }
662 public void init(FieldODEStateAndDerivative<T> state0, T t) {
663 maxError = pb.getField().getZero();
664 }
665 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
666
667 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
668 T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
669 T do1 = current.getPrimaryState()[0].subtract(theoreticalY[0]);
670 T do2 = current.getPrimaryState()[1].subtract(theoreticalY[1]);
671 T do3 = current.getPrimaryState()[2].subtract(theoreticalY[2]);
672 T error = do1.multiply(do1).add(do2.multiply(do2)).add(do3.multiply(do3));
673 if (error.subtract(maxError).getReal() > 0) {
674 maxError = error;
675 }
676 }
677 public void finish(FieldODEStateAndDerivative<T> finalState) {
678 Assert.assertEquals(0.0, maxError.getReal(), epsilon);
679 }
680 }
681
682 @Test
683
684
685
686 public abstract void testTorqueFreeMotion();
687
688 protected <T extends CalculusFieldElement<T>> void doTestTorqueFreeMotion(Field<T> field, double epsilonOmega, double epsilonQ) {
689
690 final T zero = field.getZero();
691 final T t0 = zero;
692 final T t1 = zero.newInstance(20.0);
693 final FieldVector3D<T> omegaBase = new FieldVector3D<>(zero.newInstance(5.0), zero.newInstance(0.0), zero.newInstance(4.0));
694 final FieldRotation<T> rBase = new FieldRotation<>(zero.newInstance(0.9), zero.newInstance(0.437),
695 zero.newInstance(0.0), zero.newInstance(0.0),
696 true);
697 final List<T> inertiaBase = Arrays.asList(zero.newInstance(3.0 / 8.0), zero.newInstance(1.0 / 2.0), zero.newInstance(5.0 / 8.0));
698 double minStep = 1.0e-10;
699 double maxStep = t1.subtract(t0).getReal();
700 double[] vecAbsoluteTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
701 double[] vecRelativeTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
702
703
704 Stream<TestFieldProblem8<T>> problems = Stream.empty();
705
706
707 problems = Stream.concat(problems,
708 permute(omegaBase).
709 map(omega -> new TestFieldProblem8<>(t0, t1, omega, rBase,
710 inertiaBase.get(0), FieldVector3D.getPlusI(field),
711 inertiaBase.get(1), FieldVector3D.getPlusJ(field),
712 inertiaBase.get(2), FieldVector3D.getPlusK(field))));
713
714
715 problems = Stream.concat(problems,
716 permute(rBase).
717 map(r -> new TestFieldProblem8<>(t0, t1, omegaBase, r,
718 inertiaBase.get(0), FieldVector3D.getPlusI(field),
719 inertiaBase.get(1), FieldVector3D.getPlusJ(field),
720 inertiaBase.get(2), FieldVector3D.getPlusK(field))));
721
722
723 problems = Stream.concat(problems,
724 permute(inertiaBase).
725 map(inertia -> new TestFieldProblem8<>(t0, t1, omegaBase, rBase,
726 inertia.get(0), FieldVector3D.getPlusI(field),
727 inertia.get(1), FieldVector3D.getPlusJ(field),
728 inertia.get(2), FieldVector3D.getPlusK(field))));
729
730 problems.forEach(problem -> {
731 EmbeddedRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
732 vecAbsoluteTolerance, vecRelativeTolerance);
733 integ.addStepHandler(new TorqueFreeHandler<>(problem, epsilonOmega, epsilonQ));
734 integ.integrate(new FieldExpandableODE<>(problem), problem.getInitialState(), problem.getFinalTime());
735 });
736
737 }
738
739 @Test
740 public abstract void testTorqueFreeMotionIssue230();
741
742 protected <T extends CalculusFieldElement<T>> void doTestTorqueFreeMotionIssue230(Field<T> field, double epsilonOmega, double epsilonQ) {
743
744 final T zero = field.getZero();
745 T i1 = zero.newInstance(3.0 / 8.0);
746 FieldVector3D<T> a1 = FieldVector3D.getPlusI(field);
747 T i2 = zero.newInstance(5.0 / 8.0);
748 FieldVector3D<T> a2 = FieldVector3D.getPlusK(field);
749 T i3 = zero.newInstance(1.0 / 2.0);
750 FieldVector3D<T> a3 = FieldVector3D.getPlusJ(field);
751 FieldVector3D<T> o0 = new FieldVector3D<>(zero.newInstance(5.0), zero.newInstance(0.0), zero.newInstance(4.0));
752 T o1 = FieldVector3D.dotProduct(o0, a1);
753 T o2 = FieldVector3D.dotProduct(o0, a2);
754 T o3 = FieldVector3D.dotProduct(o0, a3);
755 T e = i1.multiply(o1).multiply(o1).add(i2.multiply(o2).multiply(o2)).add(i3.multiply(o3).multiply(o3)).multiply(0.5);
756 T r1 = FastMath.sqrt(e.multiply(i1).multiply(2));
757 T r3 = FastMath.sqrt(e.multiply(i3).multiply(2));
758 int n = 50;
759 for (int i = 0; i < n; ++i) {
760 SinCos sc = FastMath.sinCos(-0.5 * FastMath.PI * (i + 50) / 200);
761 FieldVector3D<T> om = new FieldVector3D<>(r1.multiply(sc.cos()).divide(i1), a1,
762 r3.multiply(sc.sin()).divide(i3), a3);
763 TestFieldProblem8<T> problem = new TestFieldProblem8<>(zero.newInstance(0), zero.newInstance(20),
764 om,
765 new FieldRotation<>(zero.newInstance(0.9),
766 zero.newInstance(0.437),
767 zero.newInstance(0.0),
768 zero.newInstance(0.0), true),
769 i1, a1,
770 i2, a2,
771 i3, a3);
772
773 double minStep = 1.0e-10;
774 double maxStep = problem.getFinalTime().subtract(problem.getInitialTime()).getReal();
775 double[] vecAbsoluteTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
776 double[] vecRelativeTolerance = { 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14, 1.0e-14 };
777
778 EmbeddedRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
779 integ.addStepHandler(new TorqueFreeHandler<>(problem, epsilonOmega, epsilonQ));
780 integ.integrate(new FieldExpandableODE<>(problem), problem.getInitialState(), problem.getFinalTime().multiply(0.1));
781
782 }
783
784 }
785
786
787
788
789
790 private <T extends CalculusFieldElement<T>> Stream<FieldVector3D<T>> permute(final FieldVector3D<T> v) {
791 return CombinatoricsUtils.
792 permutations(Arrays.asList(v.getX(), v.getY(), v.getZ())).
793 map(a -> new FieldVector3D<>(a.get(0), a.get(1), a.get(2)));
794 }
795
796
797
798
799
800 private <T extends CalculusFieldElement<T>> Stream<FieldRotation<T>> permute(final FieldRotation<T> r) {
801 return CombinatoricsUtils.
802 permutations(Arrays.asList(r.getQ0(), r.getQ1(), r.getQ2(), r.getQ3())).
803 map(a -> new FieldRotation<>(a.get(0), a.get(1), a.get(2), a.get(3), false));
804 }
805
806
807
808
809
810 private <T extends CalculusFieldElement<T>> Stream<List<T>> permute(final List<T> list) {
811 return CombinatoricsUtils.permutations(list);
812 }
813
814 private static class TorqueFreeHandler<T extends CalculusFieldElement<T>> implements FieldODEStepHandler<T> {
815 private T maxErrorOmega;
816 private T maxErrorQ;
817 private final TestFieldProblem8<T> pb;
818 private final double epsilonOmega;
819 private final double epsilonQ;
820 private double outputStep;
821 private T current;
822
823 public TorqueFreeHandler(TestFieldProblem8<T> pb, double epsilonOmega, double epsilonQ) {
824 this.pb = pb;
825 this.epsilonOmega = epsilonOmega;
826 this.epsilonQ = epsilonQ;
827 maxErrorOmega = pb.getField().getZero();
828 maxErrorQ = pb.getField().getZero();
829 outputStep = 0.01;
830 }
831
832 public void init(FieldODEStateAndDerivative<T> state0, T t) {
833 maxErrorOmega = pb.getField().getZero();
834 maxErrorQ = pb.getField().getZero();
835 current = state0.getTime().subtract(outputStep);
836 }
837
838 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
839
840 current = current.add(outputStep);
841 while (interpolator.getPreviousState().getTime().subtract(current).getReal() <= 0 &&
842 interpolator.getCurrentState().getTime().subtract(current).getReal() > 0) {
843 FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(current);
844 final T[] theoretical = pb.computeTheoreticalState(state.getTime());
845 final T errorOmega = FieldVector3D.distance(new FieldVector3D<>(state.getPrimaryState()[0],
846 state.getPrimaryState()[1],
847 state.getPrimaryState()[2]),
848 new FieldVector3D<>(theoretical[0],
849 theoretical[1],
850 theoretical[2]));
851 maxErrorOmega = FastMath.max(maxErrorOmega, errorOmega);
852 final T errorQ = FieldRotation.distance(new FieldRotation<>(state.getPrimaryState()[3],
853 state.getPrimaryState()[4],
854 state.getPrimaryState()[5],
855 state.getPrimaryState()[6],
856 true),
857 new FieldRotation<>(theoretical[3],
858 theoretical[4],
859 theoretical[5],
860 theoretical[6],
861 true));
862 maxErrorQ = FastMath.max(maxErrorQ, errorQ);
863 current = current.add(outputStep);
864
865 }
866 }
867
868 public void finish(FieldODEStateAndDerivative<T> finalState) {
869 Assert.assertEquals(0.0, maxErrorOmega.getReal(), epsilonOmega);
870 Assert.assertEquals(0.0, maxErrorQ.getReal(), epsilonQ);
871 }
872
873 }
874
875 @Test
876 public abstract void testSecondaryEquations();
877
878 protected <T extends CalculusFieldElement<T>> void doTestSecondaryEquations(final Field<T> field,
879 final double epsilonSinCos,
880 final double epsilonLinear) {
881 FieldOrdinaryDifferentialEquation<T> sinCos = new FieldOrdinaryDifferentialEquation<T>() {
882
883 @Override
884 public int getDimension() {
885 return 2;
886 }
887
888 @Override
889 public T[] computeDerivatives(T t, T[] y) {
890 T[] yDot = y.clone();
891 yDot[0] = y[1];
892 yDot[1] = y[0].negate();
893 return yDot;
894 }
895
896 };
897
898 FieldSecondaryODE<T> linear = new FieldSecondaryODE<T>() {
899
900 @Override
901 public int getDimension() {
902 return 1;
903 }
904
905 @Override
906 public T[] computeDerivatives(T t, T[] primary, T[] primaryDot, T[] secondary) {
907 T[] secondaryDot = secondary.clone();
908 secondaryDot[0] = t.getField().getOne().negate();
909 return secondaryDot;
910 }
911
912 };
913
914 FieldExpandableODE<T> expandable = new FieldExpandableODE<>(sinCos);
915 expandable.addSecondaryEquations(linear);
916
917 FieldODEIntegrator<T> integrator = createIntegrator(field, 0.001, 1.0, 1.0e-12, 1.0e-12);
918 final double[] max = new double[2];
919 integrator.addStepHandler(new FieldODEStepHandler<T>() {
920 @Override
921 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
922 for (int i = 0; i <= 10; ++i) {
923 T tPrev = interpolator.getPreviousState().getTime();
924 T tCurr = interpolator.getCurrentState().getTime();
925 T t = tPrev.multiply(10 - i).add(tCurr.multiply(i)).divide(10);
926 FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
927 Assert.assertEquals(2, state.getPrimaryStateDimension());
928 Assert.assertEquals(1, state.getNumberOfSecondaryStates());
929 Assert.assertEquals(2, state.getSecondaryStateDimension(0));
930 Assert.assertEquals(1, state.getSecondaryStateDimension(1));
931 Assert.assertEquals(3, state.getCompleteStateDimension());
932 max[0] = FastMath.max(max[0],
933 t.sin().subtract(state.getPrimaryState()[0]).norm());
934 max[0] = FastMath.max(max[0],
935 t.cos().subtract(state.getPrimaryState()[1]).norm());
936 max[1] = FastMath.max(max[1],
937 field.getOne().subtract(t).subtract(state.getSecondaryState(1)[0]).norm());
938 }
939 }
940 });
941
942 T[] primary0 = MathArrays.buildArray(field, 2);
943 primary0[0] = field.getZero();
944 primary0[1] = field.getOne();
945 T[][] secondary0 = MathArrays.buildArray(field, 1, 1);
946 secondary0[0][0] = field.getOne();
947 FieldODEState<T> initialState = new FieldODEState<T>(field.getZero(), primary0, secondary0);
948
949 FieldODEStateAndDerivative<T> finalState =
950 integrator.integrate(expandable, initialState, field.getZero().add(10.0));
951 Assert.assertEquals(10.0, finalState.getTime().getReal(), 1.0e-12);
952 Assert.assertEquals(0, max[0], epsilonSinCos);
953 Assert.assertEquals(0, max[1], epsilonLinear);
954
955 }
956
957 @Test
958 public abstract void testPartialDerivatives();
959
960 protected void doTestPartialDerivatives(final double epsilonY,
961 final double[] epsilonPartials) {
962
963
964 final DSFactory factory = new DSFactory(5, 1);
965 final int parOmega = 0;
966 final int parTO = 1;
967 final int parY00 = 2;
968 final int parY01 = 3;
969 final int parT = 4;
970
971 DerivativeStructure omega = factory.variable(parOmega, 1.3);
972 DerivativeStructure t0 = factory.variable(parTO, 1.3);
973 DerivativeStructure[] y0 = new DerivativeStructure[] {
974 factory.variable(parY00, 3.0),
975 factory.variable(parY01, 4.0)
976 };
977 DerivativeStructure t = factory.variable(parT, 6.0);
978 SinCosODE sinCos = new SinCosODE(omega);
979
980 EmbeddedRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
981 createIntegrator(omega.getField(),
982 t.subtract(t0).multiply(0.001).getReal(), t.subtract(t0).getReal(),
983 1.0e-12, 1.0e-12);
984 FieldODEStateAndDerivative<DerivativeStructure> result =
985 integrator.integrate(new FieldExpandableODE<DerivativeStructure>(sinCos),
986 new FieldODEState<DerivativeStructure>(t0, y0),
987 t);
988
989
990 for (int i = 0; i < sinCos.getDimension(); ++i) {
991 Assert.assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getPrimaryState()[i].getValue(), epsilonY);
992 }
993
994
995 final double[][] derivatives = sinCos.getDerivatives(t.getReal());
996 for (int i = 0; i < sinCos.getDimension(); ++i) {
997 for (int parameter = 0; parameter < factory.getCompiler().getFreeParameters(); ++parameter) {
998 Assert.assertEquals(derivatives[i][parameter], dYdP(result.getPrimaryState()[i], parameter), epsilonPartials[parameter]);
999 }
1000 }
1001
1002 }
1003
1004 @Test
1005 public void testIssue118() {
1006
1007
1008 final DSFactory factory = new DSFactory(3, 3);
1009
1010
1011 final double a = 2.0;
1012 final double b = 1.0;
1013 final double omega = 0.5;
1014 final Ellipse<DerivativeStructure> ellipse =
1015 new Ellipse<>(factory.variable(0, a), factory.variable(1, b), factory.variable(2, omega));
1016 final DerivativeStructure[] initState = ellipse.computeTheoreticalState(factory.constant(0.0));
1017
1018
1019 final DerivativeStructure t0 = factory.constant(0.0);
1020 final DerivativeStructure tf = factory.constant(2.0 * FastMath.PI / omega);
1021
1022
1023 final FieldExpandableODE<DerivativeStructure> ode = new FieldExpandableODE<>(ellipse);
1024 EmbeddedRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
1025 createIntegrator(factory.getDerivativeField(), 1e-3, 1e3, 1e-12, 1e-12);
1026
1027 integrator.addStepHandler((interpolator) -> {
1028 DerivativeStructure tK = interpolator.getCurrentState().getTime();
1029 DerivativeStructure[] integrated = interpolator.getCurrentState().getPrimaryState();
1030 DerivativeStructure[] thK = ellipse.computeTheoreticalState(tK);
1031 DerivativeStructure[] tkKtrunc = ellipse.computeTheoreticalState(factory.constant(tK.getReal()));
1032 for (int i = 0 ; i < integrated.length; ++i) {
1033 final double[] integratedI = integrated[i].getAllDerivatives();
1034 final double[] theoreticalI = thK[i].getAllDerivatives();
1035 final double[] truncatedI = tkKtrunc[i].getAllDerivatives();
1036 for (int k = 0; k < factory.getCompiler().getSize(); ++k) {
1037 final int[] orders = factory.getCompiler().getPartialDerivativeOrders(k);
1038 double scaler = 1.0;
1039 for (int ord : orders) {
1040 scaler *= CombinatoricsUtils.factorialDouble(ord);
1041 }
1042 Assert.assertEquals(truncatedI[k], theoreticalI[k], 1e-15 * scaler);
1043 Assert.assertEquals(truncatedI[k], integratedI[k], 1e-8 * scaler);
1044 }
1045 }
1046 });
1047
1048 integrator.integrate(ode, new FieldODEState<>(t0, initState), tf);
1049
1050 }
1051
1052 @Test
1053 public void testInfiniteIntegration() {
1054 Field<Binary64> field = Binary64Field.getInstance();
1055 EmbeddedRungeKuttaFieldIntegrator<Binary64> fieldIntegrator = createIntegrator(Binary64Field.getInstance(), 0.01, 1.0, 0.1, 0.1);
1056 TestFieldProblem1<Binary64> pb = new TestFieldProblem1<Binary64>(field);
1057 double convergence = 1e-6;
1058 fieldIntegrator.addEventDetector(new FieldODEEventDetector<Binary64>() {
1059 @Override
1060 public FieldAdaptableInterval<Binary64> getMaxCheckInterval() {
1061 return s -> Double.POSITIVE_INFINITY;
1062 }
1063 @Override
1064 public int getMaxIterationCount() {
1065 return 1000;
1066 }
1067 @Override
1068 public BracketedRealFieldUnivariateSolver<Binary64> getSolver() {
1069 return new FieldBracketingNthOrderBrentSolver<>(new Binary64(0),
1070 new Binary64(convergence),
1071 new Binary64(0),
1072 5);
1073 }
1074 @Override
1075 public Binary64 g(FieldODEStateAndDerivative<Binary64> state) {
1076 return state.getTime().subtract(pb.getFinalTime());
1077 }
1078 @Override
1079 public FieldODEEventHandler<Binary64> getHandler() {
1080 return (state, detector, increasing) -> Action.STOP;
1081 }
1082 });
1083 FieldODEStateAndDerivative<Binary64> finalState = fieldIntegrator.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), Binary64.POSITIVE_INFINITY);
1084 Assert.assertEquals(pb.getFinalTime().getReal(), finalState.getTime().getReal(), convergence);
1085 }
1086
1087 private double dYdP(final DerivativeStructure y, final int parameter) {
1088 int[] orders = new int[y.getFreeParameters()];
1089 orders[parameter] = 1;
1090 return y.getPartialDerivative(orders);
1091 }
1092
1093 private static class SinCosODE implements FieldOrdinaryDifferentialEquation<DerivativeStructure> {
1094
1095 private final DerivativeStructure omega;
1096 private DerivativeStructure r;
1097 private DerivativeStructure alpha;
1098
1099 private double dRdY00;
1100 private double dRdY01;
1101 private double dAlphadOmega;
1102 private double dAlphadT0;
1103 private double dAlphadY00;
1104 private double dAlphadY01;
1105
1106 protected SinCosODE(final DerivativeStructure omega) {
1107 this.omega = omega;
1108 }
1109
1110 public int getDimension() {
1111 return 2;
1112 }
1113
1114 public void init(final DerivativeStructure t0, final DerivativeStructure[] y0,
1115 final DerivativeStructure finalTime) {
1116
1117
1118
1119 final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1]));
1120
1121 this.r = r2.sqrt();
1122 this.dRdY00 = y0[0].divide(r).getReal();
1123 this.dRdY01 = y0[1].divide(r).getReal();
1124
1125 this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega));
1126 this.dAlphadOmega = -t0.getReal();
1127 this.dAlphadT0 = -omega.getReal();
1128 this.dAlphadY00 = y0[1].divide(r2).getReal();
1129 this.dAlphadY01 = y0[0].negate().divide(r2).getReal();
1130
1131 }
1132
1133 public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) {
1134 return new DerivativeStructure[] {
1135 omega.multiply(y[1]),
1136 omega.multiply(y[0]).negate()
1137 };
1138 }
1139
1140 public double[] theoreticalY(final double t) {
1141 final double theta = omega.getReal() * t + alpha.getReal();
1142 return new double[] {
1143 r.getReal() * FastMath.sin(theta), r.getReal() * FastMath.cos(theta)
1144 };
1145 }
1146
1147 public double[][] getDerivatives(final double t) {
1148
1149
1150 final double theta = omega.getReal() * t + alpha.getReal();
1151 final double sin = FastMath.sin(theta);
1152 final double cos = FastMath.cos(theta);
1153 final double y0 = r.getReal() * sin;
1154 final double y1 = r.getReal() * cos;
1155
1156
1157 final double dY0dOmega = y1 * (t + dAlphadOmega);
1158 final double dY0dT0 = y1 * dAlphadT0;
1159 final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00;
1160 final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01;
1161 final double dY0dT = y1 * omega.getReal();
1162
1163
1164 final double dY1dOmega = - y0 * (t + dAlphadOmega);
1165 final double dY1dT0 = - y0 * dAlphadT0;
1166 final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00;
1167 final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01;
1168 final double dY1dT = - y0 * omega.getReal();
1169
1170 return new double[][] {
1171 { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT },
1172 { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT }
1173 };
1174
1175 }
1176
1177 }
1178
1179 }