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