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 org.hipparchus.CalculusFieldElement;
22 import org.hipparchus.Field;
23 import org.hipparchus.analysis.differentiation.DSFactory;
24 import org.hipparchus.analysis.differentiation.DerivativeStructure;
25 import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
26 import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.ode.AbstractFieldIntegrator;
30 import org.hipparchus.ode.FieldExpandableODE;
31 import org.hipparchus.ode.FieldODEIntegrator;
32 import org.hipparchus.ode.FieldODEState;
33 import org.hipparchus.ode.FieldODEStateAndDerivative;
34 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
35 import org.hipparchus.ode.FieldSecondaryODE;
36 import org.hipparchus.ode.LocalizedODEFormats;
37 import org.hipparchus.ode.MultistepFieldIntegrator;
38 import org.hipparchus.ode.TestFieldProblem1;
39 import org.hipparchus.ode.TestFieldProblem5;
40 import org.hipparchus.ode.TestFieldProblem6;
41 import org.hipparchus.ode.TestFieldProblemAbstract;
42 import org.hipparchus.ode.TestFieldProblemHandler;
43 import org.hipparchus.ode.events.Action;
44 import org.hipparchus.ode.events.FieldAdaptableInterval;
45 import org.hipparchus.ode.events.FieldODEEventDetector;
46 import org.hipparchus.ode.events.FieldODEEventHandler;
47 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
48 import org.hipparchus.ode.sampling.FieldODEStepHandler;
49 import org.hipparchus.util.Binary64Field;
50 import org.hipparchus.util.FastMath;
51 import org.hipparchus.util.MathArrays;
52 import org.junit.Assert;
53 import org.junit.Test;
54
55 public abstract class AdamsFieldIntegratorAbstractTest {
56
57 protected abstract <T extends CalculusFieldElement<T>> AdamsFieldIntegrator<T>
58 createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
59 final double scalAbsoluteTolerance, final double scalRelativeTolerance);
60
61 protected abstract <T extends CalculusFieldElement<T>> AdamsFieldIntegrator<T>
62 createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
63 final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
64
65 @Test(expected=MathIllegalArgumentException.class)
66 public abstract void testMinStep();
67
68 protected <T extends CalculusFieldElement<T>> void doDimensionCheck(final Field<T> field) {
69 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
70
71 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
72 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
73 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
74 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
75
76 FieldODEIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
77 vecAbsoluteTolerance,
78 vecRelativeTolerance);
79 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
80 integ.addStepHandler(handler);
81 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
82
83 }
84
85 protected void doNbPointsTest() {
86 try {
87 createIntegrator(Binary64Field.getInstance(), 1, 1.0e-3, 1.0e+3, 1.0e-15, 1.0e-15);
88 Assert.fail("an exception should have been thrown");
89 } catch (MathIllegalArgumentException miae) {
90 Assert.assertEquals(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
91 miae.getSpecifier());
92 }
93 try {
94 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
95 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
96 createIntegrator(Binary64Field.getInstance(),
97 1, 1.0e-3, 1.0e+3, vecAbsoluteTolerance, vecRelativeTolerance);
98 Assert.fail("an exception should have been thrown");
99 } catch (MathIllegalArgumentException miae) {
100 Assert.assertEquals(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
101 miae.getSpecifier());
102 }
103 }
104
105 @Test
106 public abstract void testIncreasingTolerance();
107
108 protected <T extends CalculusFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
109 double ratioMin, double ratioMax) {
110
111 int previousCalls = Integer.MAX_VALUE;
112 for (int i = -12; i < -2; ++i) {
113 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
114 double minStep = 0;
115 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
116 double scalAbsoluteTolerance = FastMath.pow(10.0, i);
117 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
118
119 MultistepFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
120 scalAbsoluteTolerance,
121 scalRelativeTolerance);
122 int orderCorrection = integ instanceof AdamsBashforthFieldIntegrator ? 0 : 1;
123 Assert.assertEquals(FastMath.pow(2.0, 1.0 / (4 + orderCorrection)), integ.getMaxGrowth(), 1.0e-10);
124 Assert.assertEquals(0.2, integ.getMinReduction(), 1.0e-10);
125 Assert.assertEquals(4, integ.getNSteps());
126 Assert.assertEquals(0.9, integ.getSafety(), 1.0e-10);
127 Assert.assertTrue(integ.getStarterIntegrator() instanceof DormandPrince853FieldIntegrator);
128 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
129 integ.addStepHandler(handler);
130 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
131
132 Assert.assertTrue(handler.getMaximalValueError().getReal() > ratioMin * scalAbsoluteTolerance);
133 Assert.assertTrue(handler.getMaximalValueError().getReal() < ratioMax * scalAbsoluteTolerance);
134
135 int calls = pb.getCalls();
136 Assert.assertEquals(integ.getEvaluations(), calls);
137 Assert.assertTrue(calls <= previousCalls);
138 previousCalls = calls;
139
140 }
141
142 }
143
144 @Test(expected = MathIllegalStateException.class)
145 public abstract void exceedMaxEvaluations();
146
147 protected <T extends CalculusFieldElement<T>> void doExceedMaxEvaluations(final Field<T> field, final int max) {
148
149 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
150 double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
151
152 FieldODEIntegrator<T> integ = createIntegrator(field, 2, 0, range, 1.0e-12, 1.0e-12);
153 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
154 integ.addStepHandler(handler);
155 integ.setMaxEvaluations(max);
156 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
157
158 }
159
160 @Test
161 public abstract void backward();
162
163 protected <T extends CalculusFieldElement<T>> void doBackward(final Field<T> field,
164 final double epsilonLast,
165 final double epsilonMaxValue,
166 final double epsilonMaxTime,
167 final String name) {
168
169 final double resetTime = -3.98;
170 final TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field) {
171 @Override
172 public T[] getTheoreticalEventsTimes() {
173 final T[] tEv = MathArrays.buildArray(field, 1);
174 tEv[0] = field.getZero().add(resetTime);
175 return tEv;
176 }
177 };
178 double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
179
180 AdamsFieldIntegrator<T> integ = createIntegrator(field, 4, 0, range, 1.0e-12, 1.0e-12);
181 FieldODEEventDetector<T> event = new FieldODEEventDetector<T>() {
182
183 @Override
184 public FieldAdaptableInterval<T> getMaxCheckInterval() {
185 return s -> 0.5 * range;
186 }
187
188 @Override
189 public int getMaxIterationCount() {
190 return 100;
191 }
192
193 @Override
194 public BracketedRealFieldUnivariateSolver<T> getSolver() {
195 return new FieldBracketingNthOrderBrentSolver<T>(field.getZero(),
196 field.getZero().newInstance(1.0e-6 * range),
197 field.getZero(),
198 5);
199 }
200
201 @Override
202 public FieldODEEventHandler<T> getHandler() {
203 return (state, detector, increasing) -> Action.RESET_STATE;
204 }
205 @Override
206 public T g(FieldODEStateAndDerivative<T> state) {
207 return state.getTime().subtract(resetTime);
208 }
209
210 };
211 integ.addEventDetector(event);
212 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
213 integ.addStepHandler(handler);
214 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
215
216 Assert.assertEquals(0.0, handler.getLastError().getReal(), epsilonLast);
217 Assert.assertEquals(0.0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
218 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
219 Assert.assertEquals(name, integ.getName());
220 }
221
222 @Test
223 public abstract void polynomial();
224
225 protected <T extends CalculusFieldElement<T>> void doPolynomial(final Field<T> field,
226 final int nLimit,
227 final double epsilonBad,
228 final double epsilonGood) {
229 final TestFieldProblem6<T> pb = new TestFieldProblem6<T>(field);
230 final double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).norm();
231
232 for (int nSteps = 2; nSteps < 8; ++nSteps) {
233 AdamsFieldIntegrator<T> integ = createIntegrator(field, nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
234 FieldODEEventDetector<T> event = new FieldODEEventDetector<T>() {
235
236 @Override
237 public FieldAdaptableInterval<T> getMaxCheckInterval() {
238 return s -> 0.5 * range;
239 }
240
241 @Override
242 public int getMaxIterationCount() {
243 return 100;
244 }
245
246 @Override
247 public BracketedRealFieldUnivariateSolver<T> getSolver() {
248 return new FieldBracketingNthOrderBrentSolver<T>(field.getZero(),
249 field.getZero().newInstance(1.0e-6 * range),
250 field.getZero(),
251 5);
252 }
253
254 @Override
255 public FieldODEEventHandler<T> getHandler() {
256 return (state, detector, increasing) -> Action.RESET_STATE;
257 }
258
259
260 @Override
261 public T g(FieldODEStateAndDerivative<T> state) {
262 return state.getTime().subtract(pb.getInitialState().getTime().getReal() + 0.5 * range);
263 }
264
265 };
266 integ.addEventDetector(event);
267 integ.setStarterIntegrator(new PerfectStarter<T>(pb, nSteps));
268 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
269 integ.addStepHandler(handler);
270 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
271 if (nSteps < nLimit) {
272 Assert.assertTrue(handler.getMaximalValueError().getReal() > epsilonBad);
273 } else {
274 Assert.assertTrue(handler.getMaximalValueError().getReal() < epsilonGood);
275 }
276 }
277
278 }
279
280 @Test
281 public void testNaNAppearing() {
282 doTestNaNAppearing(Binary64Field.getInstance());
283 }
284
285 private <T extends CalculusFieldElement<T>> void doTestNaNAppearing(final Field<T> field) {
286 try {
287 AdamsFieldIntegrator<T> integ = createIntegrator(field, 8, 0.01, 1.0, 0.1, 0.1);
288 final FieldOrdinaryDifferentialEquation<T> ode = new FieldOrdinaryDifferentialEquation<T>() {
289 public int getDimension() {
290 return 1;
291 }
292 public T[] computeDerivatives(T t, T[] y) {
293 T[] yDot = MathArrays.buildArray(t.getField(), getDimension());
294 yDot[0] = FastMath.log(t);
295 return yDot;
296 }
297 };
298 final T t0 = field.getZero().add(10.0);
299 final T t1 = field.getZero().add(-1.0);
300 final T[] y0 = MathArrays.buildArray(field, ode.getDimension());
301 y0[0] = field.getZero().add(1.0);
302 integ.integrate(new FieldExpandableODE<>(ode), new FieldODEState<>(t0, y0), t1);
303 Assert.fail("an exception should have been thrown");
304 } catch (MathIllegalStateException mise) {
305 Assert.assertEquals(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION, mise.getSpecifier());
306 Assert.assertTrue(((Double) mise.getParts()[0]).doubleValue() <= 0.0);
307 }
308 }
309
310 @Test
311 public abstract void testSecondaryEquations();
312
313 protected <T extends CalculusFieldElement<T>> void doTestSecondaryEquations(final Field<T> field,
314 final double epsilonSinCos,
315 final double epsilonLinear) {
316 FieldOrdinaryDifferentialEquation<T> sinCos = new FieldOrdinaryDifferentialEquation<T>() {
317
318 @Override
319 public int getDimension() {
320 return 2;
321 }
322
323 @Override
324 public T[] computeDerivatives(T t, T[] y) {
325 T[] yDot = y.clone();
326 yDot[0] = y[1];
327 yDot[1] = y[0].negate();
328 return yDot;
329 }
330
331 };
332
333 FieldSecondaryODE<T> linear = new FieldSecondaryODE<T>() {
334
335 @Override
336 public int getDimension() {
337 return 1;
338 }
339
340 @Override
341 public T[] computeDerivatives(T t, T[] primary, T[] primaryDot, T[] secondary) {
342 T[] secondaryDot = secondary.clone();
343 secondaryDot[0] = t.getField().getOne().negate();
344 return secondaryDot;
345 }
346
347 };
348
349 FieldExpandableODE<T> expandable = new FieldExpandableODE<>(sinCos);
350 expandable.addSecondaryEquations(linear);
351
352 FieldODEIntegrator<T> integrator = createIntegrator(field, 6, 0.001, 1.0, 1.0e-12, 1.0e-12);
353 final double[] max = new double[2];
354 integrator.addStepHandler(new FieldODEStepHandler<T>() {
355 @Override
356 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
357 for (int i = 0; i <= 10; ++i) {
358 T tPrev = interpolator.getPreviousState().getTime();
359 T tCurr = interpolator.getCurrentState().getTime();
360 T t = tPrev.multiply(10 - i).add(tCurr.multiply(i)).divide(10);
361 FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
362 Assert.assertEquals(2, state.getPrimaryStateDimension());
363 Assert.assertEquals(1, state.getNumberOfSecondaryStates());
364 Assert.assertEquals(2, state.getSecondaryStateDimension(0));
365 Assert.assertEquals(1, state.getSecondaryStateDimension(1));
366 Assert.assertEquals(3, state.getCompleteStateDimension());
367 max[0] = FastMath.max(max[0],
368 t.sin().subtract(state.getPrimaryState()[0]).norm());
369 max[0] = FastMath.max(max[0],
370 t.cos().subtract(state.getPrimaryState()[1]).norm());
371 max[1] = FastMath.max(max[1],
372 field.getOne().subtract(t).subtract(state.getSecondaryState(1)[0]).norm());
373 }
374 }
375 });
376
377 T[] primary0 = MathArrays.buildArray(field, 2);
378 primary0[0] = field.getZero();
379 primary0[1] = field.getOne();
380 T[][] secondary0 = MathArrays.buildArray(field, 1, 1);
381 secondary0[0][0] = field.getOne();
382 FieldODEState<T> initialState = new FieldODEState<T>(field.getZero(), primary0, secondary0);
383
384 FieldODEStateAndDerivative<T> finalState =
385 integrator.integrate(expandable, initialState, field.getZero().add(10.0));
386 Assert.assertEquals(10.0, finalState.getTime().getReal(), 1.0e-12);
387 Assert.assertEquals(0, max[0], epsilonSinCos);
388 Assert.assertEquals(0, max[1], epsilonLinear);
389
390 }
391
392 @Test(expected=MathIllegalStateException.class)
393 public abstract void testStartFailure();
394
395 protected <T extends CalculusFieldElement<T>> void doTestStartFailure(final Field<T> field) {
396 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
397 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0001).getReal();
398 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
399 double scalAbsoluteTolerance = 1.0e-6;
400 double scalRelativeTolerance = 1.0e-7;
401
402 MultistepFieldIntegrator<T> integ = createIntegrator(field, 6, minStep, maxStep,
403 scalAbsoluteTolerance,
404 scalRelativeTolerance);
405 integ.setStarterIntegrator(new DormandPrince853FieldIntegrator<T>(field, maxStep * 0.5, maxStep, 0.1, 0.1));
406 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
407 integ.addStepHandler(handler);
408 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
409
410 }
411
412 @Test
413 public void testIssue118() {
414
415
416 final DSFactory factory = new DSFactory(3, 3);
417
418
419 final double a = 2.0;
420 final double b = 1.0;
421 final double omega = 0.5;
422 final Ellipse<DerivativeStructure> ellipse =
423 new Ellipse<>(factory.variable(0, a), factory.variable(1, b), factory.variable(2, omega));
424 final DerivativeStructure[] initState = ellipse.computeTheoreticalState(factory.constant(0.0));
425
426
427 final DerivativeStructure t0 = factory.constant(0.0);
428 final DerivativeStructure tf = factory.constant(2.0 * FastMath.PI / omega);
429
430
431 final FieldExpandableODE<DerivativeStructure> ode = new FieldExpandableODE<>(ellipse);
432 MultistepFieldIntegrator<DerivativeStructure> integrator =
433 createIntegrator(factory.getDerivativeField(), 6, 1e-3, 1e3, 1e-12, 1e-12);
434
435 integrator.addStepHandler((interpolator) -> {
436 DerivativeStructure tK = interpolator.getCurrentState().getTime();
437 DerivativeStructure[] integrated = interpolator.getCurrentState().getPrimaryState();
438 DerivativeStructure[] thK = ellipse.computeTheoreticalState(tK);
439 DerivativeStructure[] tkKtrunc = ellipse.computeTheoreticalState(factory.constant(tK.getReal()));
440 for (int i = 0 ; i < integrated.length; ++i) {
441 final double[] integratedI = integrated[i].getAllDerivatives();
442 final double[] theoreticalI = thK[i].getAllDerivatives();
443 final double[] truncatedI = tkKtrunc[i].getAllDerivatives();
444 for (int k = 0; k < factory.getCompiler().getSize(); ++k) {
445 Assert.assertEquals(truncatedI[k], theoreticalI[k], 1e-15);
446 Assert.assertEquals(truncatedI[k], integratedI[k], 3e-6);
447 }
448 }
449 });
450
451 integrator.integrate(ode, new FieldODEState<>(t0, initState), tf);
452
453 }
454
455 private static class PerfectStarter<T extends CalculusFieldElement<T>> extends AbstractFieldIntegrator<T> {
456
457 private final PerfectInterpolator<T> interpolator;
458 private final int nbSteps;
459
460 public PerfectStarter(final TestFieldProblemAbstract<T> problem, final int nbSteps) {
461 super(problem.getField(), "perfect-starter");
462 this.interpolator = new PerfectInterpolator<T>(problem);
463 this.nbSteps = nbSteps;
464 }
465
466 public FieldODEStateAndDerivative<T> integrate(FieldExpandableODE<T> equations,
467 FieldODEState<T> initialState, T finalTime) {
468 T tStart = initialState.getTime().add(finalTime.subtract(initialState.getTime()).multiply(0.01));
469 getEvaluationsCounter().increment(nbSteps);
470 interpolator.setCurrentTime(initialState.getTime());
471 for (int i = 0; i < nbSteps; ++i) {
472 T tK = initialState.getTime().multiply(nbSteps - 1 - (i + 1)).add(tStart.multiply(i + 1)).divide(nbSteps - 1);
473 interpolator.setPreviousTime(interpolator.getCurrentTime());
474 interpolator.setCurrentTime(tK);
475 for (FieldODEStepHandler<T> handler : getStepHandlers()) {
476 handler.handleStep(interpolator);
477 if (i == nbSteps - 1) {
478 handler.finish(interpolator.getCurrentState());
479 }
480 }
481 }
482 return interpolator.getInterpolatedState(tStart);
483 }
484
485 }
486
487 private static class PerfectInterpolator<T extends CalculusFieldElement<T>> implements FieldODEStateInterpolator<T> {
488 private final TestFieldProblemAbstract<T> problem;
489 private T previousTime;
490 private T currentTime;
491
492 public PerfectInterpolator(final TestFieldProblemAbstract<T> problem) {
493 this.problem = problem;
494 }
495
496 public void setPreviousTime(T previousTime) {
497 this.previousTime = previousTime;
498 }
499
500 public void setCurrentTime(T currentTime) {
501 this.currentTime = currentTime;
502 }
503
504 public T getCurrentTime() {
505 return currentTime;
506 }
507
508 public boolean isForward() {
509 return problem.getFinalTime().subtract(problem.getInitialState().getTime()).getReal() >= 0;
510 }
511
512 public FieldODEStateAndDerivative<T> getPreviousState() {
513 return getInterpolatedState(previousTime);
514 }
515
516 @Override
517 public boolean isPreviousStateInterpolated() {
518 return false;
519 }
520
521 public FieldODEStateAndDerivative<T> getCurrentState() {
522 return getInterpolatedState(currentTime);
523 }
524
525 @Override
526 public boolean isCurrentStateInterpolated() {
527 return false;
528 }
529
530 public FieldODEStateAndDerivative<T> getInterpolatedState(T time) {
531 T[] y = problem.computeTheoreticalState(time);
532 T[] yDot = problem.computeDerivatives(time, y);
533 return new FieldODEStateAndDerivative<T>(time, y, yDot);
534 }
535
536 }
537
538 }