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