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