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.analysis.UnivariateFunction;
22 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
23 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
24 import org.hipparchus.exception.MathIllegalArgumentException;
25 import org.hipparchus.exception.MathIllegalStateException;
26 import org.hipparchus.ode.AbstractIntegrator;
27 import org.hipparchus.ode.ExpandableODE;
28 import org.hipparchus.ode.LocalizedODEFormats;
29 import org.hipparchus.ode.MultistepIntegrator;
30 import org.hipparchus.ode.ODEIntegrator;
31 import org.hipparchus.ode.ODEState;
32 import org.hipparchus.ode.ODEStateAndDerivative;
33 import org.hipparchus.ode.OrdinaryDifferentialEquation;
34 import org.hipparchus.ode.SecondaryODE;
35 import org.hipparchus.ode.TestProblem1;
36 import org.hipparchus.ode.TestProblem5;
37 import org.hipparchus.ode.TestProblem6;
38 import org.hipparchus.ode.TestProblemAbstract;
39 import org.hipparchus.ode.TestProblemHandler;
40 import org.hipparchus.ode.events.Action;
41 import org.hipparchus.ode.events.AdaptableInterval;
42 import org.hipparchus.ode.events.ODEEventDetector;
43 import org.hipparchus.ode.events.ODEEventHandler;
44 import org.hipparchus.ode.sampling.ODEStateInterpolator;
45 import org.hipparchus.ode.sampling.ODEStepHandler;
46 import org.hipparchus.util.FastMath;
47 import org.junit.jupiter.api.Test;
48
49 import static org.junit.jupiter.api.Assertions.assertEquals;
50 import static org.junit.jupiter.api.Assertions.assertTrue;
51 import static org.junit.jupiter.api.Assertions.fail;
52
53 public abstract class AdamsIntegratorAbstractTest {
54
55 protected abstract AdamsIntegrator
56 createIntegrator(final int nSteps, final double minStep, final double maxStep,
57 final double scalAbsoluteTolerance, final double scalRelativeTolerance);
58
59 protected abstract AdamsIntegrator
60 createIntegrator(final int nSteps, final double minStep, final double maxStep,
61 final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
62
63 @Test
64 public abstract void testMinStep();
65
66 protected void doNbPointsTest() {
67 try {
68 createIntegrator(1, 1.0e-3, 1.0e+3, 1.0e-15, 1.0e-15);
69 fail("an exception should have been thrown");
70 } catch (MathIllegalArgumentException miae) {
71 assertEquals(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
72 miae.getSpecifier());
73 }
74 try {
75 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
76 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
77 createIntegrator(1, 1.0e-3, 1.0e+3, vecAbsoluteTolerance, vecRelativeTolerance);
78 fail("an exception should have been thrown");
79 } catch (MathIllegalArgumentException miae) {
80 assertEquals(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
81 miae.getSpecifier());
82 }
83 }
84
85 protected void doDimensionCheck() {
86 TestProblem1 pb = new TestProblem1();
87
88 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialState().getTime());
89 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
90 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
91 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
92
93 ODEIntegrator integ = createIntegrator(4, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
94 TestProblemHandler handler = new TestProblemHandler(pb, integ);
95 integ.addStepHandler(handler);
96 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
97
98 }
99
100 @Test
101 public abstract void testIncreasingTolerance();
102
103 protected void doTestIncreasingTolerance(double ratioMin, double ratioMax) {
104
105 int previousCalls = Integer.MAX_VALUE;
106 for (int i = -12; i < -2; ++i) {
107 TestProblem1 pb = new TestProblem1();
108 double minStep = 0;
109 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
110 double scalAbsoluteTolerance = FastMath.pow(10.0, i);
111 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
112
113 MultistepIntegrator integ = createIntegrator(4, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
114 int orderCorrection = integ instanceof AdamsBashforthIntegrator ? 0 : 1;
115 assertEquals(FastMath.pow(2.0, 1.0 / (4 + orderCorrection)), integ.getMaxGrowth(), 1.0e-10);
116 assertEquals(0.2, integ.getMinReduction(), 1.0e-10);
117 assertEquals(4, integ.getNSteps());
118 assertEquals(0.9, integ.getSafety(), 1.0e-10);
119 assertTrue(integ.getStarterIntegrator() instanceof DormandPrince853Integrator);
120 TestProblemHandler handler = new TestProblemHandler(pb, integ);
121 integ.addStepHandler(handler);
122 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
123
124 assertTrue(handler.getMaximalValueError() > ratioMin * scalAbsoluteTolerance);
125 assertTrue(handler.getMaximalValueError() < ratioMax * scalAbsoluteTolerance);
126
127 int calls = pb.getCalls();
128 assertEquals(integ.getEvaluations(), calls);
129 assertTrue(calls <= previousCalls);
130 previousCalls = calls;
131
132 }
133
134 }
135
136 @Test
137 public abstract void exceedMaxEvaluations();
138
139 protected void doExceedMaxEvaluations(final int max) {
140
141 TestProblem1 pb = new TestProblem1();
142 double range = pb.getFinalTime() - pb.getInitialState().getTime();
143
144 ODEIntegrator integ = createIntegrator(2, 0, range, 1.0e-12, 1.0e-12);
145 TestProblemHandler handler = new TestProblemHandler(pb, integ);
146 integ.addStepHandler(handler);
147 integ.setMaxEvaluations(max);
148 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
149
150 }
151
152 @Test
153 public abstract void backward();
154
155 protected void doBackward(final double epsilonLast,
156 final double epsilonMaxValue,
157 final double epsilonMaxTime,
158 final String name) {
159
160 final double resetTime = -3.98;
161 final TestProblem5 pb = new TestProblem5() {
162 @Override
163 public double[] getTheoreticalEventsTimes() {
164 return new double[] { resetTime };
165 }
166 };
167 final double range = pb.getFinalTime() - pb.getInitialState().getTime();
168
169 AdamsIntegrator integ = createIntegrator(4, 0, range, 1.0e-12, 1.0e-12);
170 ODEEventDetector event = new ODEEventDetector() {
171
172 @Override
173 public AdaptableInterval getMaxCheckInterval() {
174 return (s, isForward) -> 0.5 * range;
175 }
176
177 @Override
178 public int getMaxIterationCount() {
179 return 100;
180 }
181
182 @Override
183 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
184 return new BracketingNthOrderBrentSolver(0, 1.0e-6 * range, 0, 5);
185 }
186
187 @Override
188 public ODEEventHandler getHandler() {
189 return (state, detector, increasing) -> Action.RESET_STATE;
190 }
191
192 @Override
193 public double g(ODEStateAndDerivative state) {
194 return state.getTime() - resetTime;
195 }
196
197 };
198 integ.addEventDetector(event);
199 TestProblemHandler handler = new TestProblemHandler(pb, integ);
200 integ.addStepHandler(handler);
201 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
202
203 assertEquals(0.0, handler.getLastError(), epsilonLast);
204 assertEquals(0.0, handler.getMaximalValueError(), epsilonMaxValue);
205 assertEquals(0, handler.getMaximalTimeError(), epsilonMaxTime);
206 assertEquals(name, integ.getName());
207 }
208
209 @Test
210 public abstract void polynomial();
211
212 protected void doPolynomial(final int nLimit, final double epsilonBad, final double epsilonGood) {
213 TestProblem6 pb = new TestProblem6();
214 double range = FastMath.abs(pb.getFinalTime() - pb.getInitialState().getTime());
215
216 for (int nSteps = 2; nSteps < 8; ++nSteps) {
217 AdamsIntegrator integ = createIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
218 ODEEventDetector event = new ODEEventDetector() {
219
220 @Override
221 public AdaptableInterval getMaxCheckInterval() {
222 return (s, isForward) -> 0.5 * range;
223 }
224
225 @Override
226 public int getMaxIterationCount() {
227 return 100;
228 }
229
230 @Override
231 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
232 return new BracketingNthOrderBrentSolver(0, 1.0e-6 * range, 0, 5);
233 }
234
235 @Override
236 public ODEEventHandler getHandler() {
237 return (state, detector, increasing) -> Action.RESET_STATE;
238 }
239
240 @Override
241 public double g(ODEStateAndDerivative state) {
242 return state.getTime() - (pb.getInitialState().getTime() + 0.5 * range);
243 }
244
245 };
246 integ.addEventDetector(event);
247 integ.setStarterIntegrator(new PerfectStarter(pb, nSteps));
248 TestProblemHandler handler = new TestProblemHandler(pb, integ);
249 integ.addStepHandler(handler);
250 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
251 if (nSteps < nLimit) {
252 assertTrue(handler.getMaximalValueError() > epsilonBad);
253 } else {
254 assertTrue(handler.getMaximalValueError() < epsilonGood);
255 }
256 }
257
258 }
259
260 @Test
261 public void testNaNAppearing() {
262 try {
263 ODEIntegrator integ = createIntegrator(8, 0.01, 1.0, 0.1, 0.1);
264 integ.integrate(new OrdinaryDifferentialEquation() {
265 public int getDimension() {
266 return 1;
267 }
268 public double[] computeDerivatives(double t, double[] y) {
269 return new double[] { FastMath.log(t) };
270 }
271 }, new ODEState(10.0, new double[] { 1.0 }), -1.0);
272 fail("an exception should have been thrown");
273 } catch (MathIllegalStateException mise) {
274 assertEquals(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION, mise.getSpecifier());
275 assertTrue(((Double) mise.getParts()[0]).doubleValue() <= 0.0);
276 }
277 }
278
279 @Test
280 public abstract void testSecondaryEquations();
281
282 protected void doTestSecondaryEquations(final double epsilonSinCos,
283 final double epsilonLinear) {
284 OrdinaryDifferentialEquation sinCos = new OrdinaryDifferentialEquation() {
285
286 @Override
287 public int getDimension() {
288 return 2;
289 }
290
291 @Override
292 public double[] computeDerivatives(double t, double[] y) {
293 return new double[] { y[1], -y[0] };
294 }
295
296 };
297
298 SecondaryODE linear = new SecondaryODE() {
299
300 @Override
301 public int getDimension() {
302 return 1;
303 }
304
305 @Override
306 public double[] computeDerivatives(double t, double[] primary, double[] primaryDot, double[] secondary) {
307 return new double[] { -1 };
308 }
309
310 };
311
312 ExpandableODE expandable = new ExpandableODE(sinCos);
313 expandable.addSecondaryEquations(linear);
314
315 ODEIntegrator integrator = createIntegrator(6, 0.001, 1.0, 1.0e-12, 1.0e-12);
316 final double[] max = new double[2];
317 integrator.addStepHandler(new ODEStepHandler() {
318 @Override
319 public void handleStep(ODEStateInterpolator interpolator) {
320 for (int i = 0; i <= 10; ++i) {
321 double tPrev = interpolator.getPreviousState().getTime();
322 double tCurr = interpolator.getCurrentState().getTime();
323 double t = (tPrev * (10 - i) + tCurr * i) / 10;
324 ODEStateAndDerivative state = interpolator.getInterpolatedState(t);
325 assertEquals(2, state.getPrimaryStateDimension());
326 assertEquals(1, state.getNumberOfSecondaryStates());
327 assertEquals(2, state.getSecondaryStateDimension(0));
328 assertEquals(1, state.getSecondaryStateDimension(1));
329 assertEquals(3, state.getCompleteStateDimension());
330 max[0] = FastMath.max(max[0],
331 FastMath.abs(FastMath.sin(t) - state.getPrimaryState()[0]));
332 max[0] = FastMath.max(max[0],
333 FastMath.abs(FastMath.cos(t) - state.getPrimaryState()[1]));
334 max[1] = FastMath.max(max[1],
335 FastMath.abs(1 - t - state.getSecondaryState(1)[0]));
336 }
337 }
338 });
339
340 double[] primary0 = new double[] { 0.0, 1.0 };
341 double[][] secondary0 = new double[][] { { 1.0 } };
342 ODEState initialState = new ODEState(0.0, primary0, secondary0);
343
344 ODEStateAndDerivative finalState =
345 integrator.integrate(expandable, initialState, 10.0);
346 assertEquals(10.0, finalState.getTime(), 1.0e-12);
347 assertEquals(0, max[0], epsilonSinCos);
348 assertEquals(0, max[1], epsilonLinear);
349
350 }
351
352 @Test
353 public abstract void testStartFailure();
354
355 protected void doTestStartFailure() {
356 TestProblem1 pb = new TestProblem1();
357 double minStep = 0.0001 * (pb.getFinalTime() - pb.getInitialState().getTime());
358 double maxStep = pb.getFinalTime() - pb.getInitialState().getTime();
359 double scalAbsoluteTolerance = 1.0e-6;
360 double scalRelativeTolerance = 1.0e-7;
361
362 MultistepIntegrator integ = createIntegrator(6, minStep, maxStep,
363 scalAbsoluteTolerance, scalRelativeTolerance);
364 integ.setStarterIntegrator(new DormandPrince853Integrator(maxStep * 0.5, maxStep, 0.1, 0.1));
365 TestProblemHandler handler = new TestProblemHandler(pb, integ);
366 integ.addStepHandler(handler);
367 integ.integrate(new ExpandableODE(pb), pb.getInitialState(), pb.getFinalTime());
368
369 }
370
371 private static class PerfectStarter extends AbstractIntegrator {
372
373 private final PerfectInterpolator interpolator;
374 private final int nbSteps;
375
376 public PerfectStarter(final TestProblemAbstract problem, final int nbSteps) {
377 super("perfect-starter");
378 this.interpolator = new PerfectInterpolator(problem);
379 this.nbSteps = nbSteps;
380 }
381
382 public ODEStateAndDerivative integrate(ExpandableODE equations,
383 ODEState initialState, double finalTime) {
384 double tStart = initialState.getTime() + 0.01 * (finalTime - initialState.getTime());
385 getEvaluationsCounter().increment(nbSteps);
386 interpolator.setCurrentTime(initialState.getTime());
387 for (int i = 0; i < nbSteps; ++i) {
388 double tK = ((nbSteps - 1 - (i + 1)) * initialState.getTime() + (i + 1) * tStart) /
389 (nbSteps - 1);
390 interpolator.setPreviousTime(interpolator.getCurrentTime());
391 interpolator.setCurrentTime(tK);
392 for (ODEStepHandler handler : getStepHandlers()) {
393 handler.handleStep(interpolator);
394 if (i == nbSteps - 1) {
395 handler.finish(interpolator.getCurrentState());
396 }
397 }
398 }
399 return interpolator.getInterpolatedState(tStart);
400 }
401
402 }
403
404 private static class PerfectInterpolator implements ODEStateInterpolator {
405 private static final long serialVersionUID = 20160417L;
406 private final TestProblemAbstract problem;
407 private double previousTime;
408 private double currentTime;
409
410 public PerfectInterpolator(final TestProblemAbstract problem) {
411 this.problem = problem;
412 }
413
414 public void setPreviousTime(double previousTime) {
415 this.previousTime = previousTime;
416 }
417
418 public void setCurrentTime(double currentTime) {
419 this.currentTime = currentTime;
420 }
421
422 public double getCurrentTime() {
423 return currentTime;
424 }
425
426 public boolean isForward() {
427 return problem.getFinalTime() >= problem.getInitialState().getTime();
428 }
429
430 public ODEStateAndDerivative getPreviousState() {
431 return getInterpolatedState(previousTime);
432 }
433
434 @Override
435 public boolean isPreviousStateInterpolated() {
436 return false;
437 }
438
439 public ODEStateAndDerivative getCurrentState() {
440 return getInterpolatedState(currentTime);
441 }
442
443 @Override
444 public boolean isCurrentStateInterpolated() {
445 return false;
446 }
447
448 public ODEStateAndDerivative getInterpolatedState(double time) {
449 double[] y = problem.computeTheoreticalState(time);
450 double[] yDot = problem.computeDerivatives(time, y);
451 return new ODEStateAndDerivative(time, y, yDot);
452 }
453
454 @Override
455 public ODEStateInterpolator restrictStep(ODEStateAndDerivative previousState, ODEStateAndDerivative currentState) {
456 return this;
457 }
458 }
459
460 }