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