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