1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20 import org.apache.commons.math.ode.DerivativeException;
21 import org.apache.commons.math.ode.FirstOrderIntegrator;
22 import org.apache.commons.math.ode.IntegratorException;
23 import org.apache.commons.math.ode.events.EventHandler;
24 import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
25 import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
26 import org.apache.commons.math.ode.sampling.StepHandler;
27 import org.apache.commons.math.ode.sampling.StepInterpolator;
28
29 import junit.framework.*;
30
31 public class DormandPrince54IntegratorTest
32 extends TestCase {
33
34 public DormandPrince54IntegratorTest(String name) {
35 super(name);
36 }
37
38 public void testDimensionCheck() {
39 try {
40 TestProblem1 pb = new TestProblem1();
41 DormandPrince54Integrator integrator = new DormandPrince54Integrator(0.0, 1.0,
42 1.0e-10, 1.0e-10);
43 integrator.integrate(pb,
44 0.0, new double[pb.getDimension()+10],
45 1.0, new double[pb.getDimension()+10]);
46 fail("an exception should have been thrown");
47 } catch(DerivativeException de) {
48 fail("wrong exception caught");
49 } catch(IntegratorException ie) {
50 }
51 }
52
53 public void testMinStep()
54 throws DerivativeException, IntegratorException {
55
56 try {
57 TestProblem1 pb = new TestProblem1();
58 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
59 double maxStep = pb.getFinalTime() - pb.getInitialTime();
60 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
61 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
62
63 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
64 vecAbsoluteTolerance,
65 vecRelativeTolerance);
66 TestProblemHandler handler = new TestProblemHandler(pb, integ);
67 integ.addStepHandler(handler);
68 integ.integrate(pb,
69 pb.getInitialTime(), pb.getInitialState(),
70 pb.getFinalTime(), new double[pb.getDimension()]);
71 fail("an exception should have been thrown");
72 } catch(DerivativeException de) {
73 fail("wrong exception caught");
74 } catch(IntegratorException ie) {
75 }
76
77 }
78
79 public void testSmallLastStep()
80 throws DerivativeException, IntegratorException {
81
82 TestProblemAbstract pb = new TestProblem5();
83 double minStep = 1.25;
84 double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
85 double scalAbsoluteTolerance = 6.0e-4;
86 double scalRelativeTolerance = 6.0e-4;
87
88 AdaptiveStepsizeIntegrator integ =
89 new DormandPrince54Integrator(minStep, maxStep,
90 scalAbsoluteTolerance,
91 scalRelativeTolerance);
92
93 DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
94 integ.addStepHandler(handler);
95 integ.setInitialStepSize(1.7);
96 integ.integrate(pb,
97 pb.getInitialTime(), pb.getInitialState(),
98 pb.getFinalTime(), new double[pb.getDimension()]);
99 assertTrue(handler.wasLastSeen());
100 assertEquals("Dormand-Prince 5(4)", integ.getName());
101
102 }
103
104 public void testBackward()
105 throws DerivativeException, IntegratorException {
106
107 TestProblem5 pb = new TestProblem5();
108 double minStep = 0;
109 double maxStep = pb.getFinalTime() - pb.getInitialTime();
110 double scalAbsoluteTolerance = 1.0e-8;
111 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
112
113 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
114 scalAbsoluteTolerance,
115 scalRelativeTolerance);
116 TestProblemHandler handler = new TestProblemHandler(pb, integ);
117 integ.addStepHandler(handler);
118 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
119 pb.getFinalTime(), new double[pb.getDimension()]);
120
121 assertTrue(handler.getLastError() < 2.0e-7);
122 assertTrue(handler.getMaximalValueError() < 2.0e-7);
123 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
124 assertEquals("Dormand-Prince 5(4)", integ.getName());
125 }
126
127 private static class DP54SmallLastHandler implements StepHandler {
128
129 private static final long serialVersionUID = -8168590945325629799L;
130
131 public DP54SmallLastHandler(double minStep) {
132 lastSeen = false;
133 this.minStep = minStep;
134 }
135
136 public boolean requiresDenseOutput() {
137 return false;
138 }
139
140 public void reset() {
141 }
142
143 public void handleStep(StepInterpolator interpolator, boolean isLast) {
144 if (isLast) {
145 lastSeen = true;
146 double h = interpolator.getCurrentTime() - interpolator.getPreviousTime();
147 assertTrue(Math.abs(h) < minStep);
148 }
149 }
150
151 public boolean wasLastSeen() {
152 return lastSeen;
153 }
154
155 private boolean lastSeen;
156 private double minStep;
157
158 }
159
160 public void testIncreasingTolerance()
161 throws DerivativeException, IntegratorException {
162
163 int previousCalls = Integer.MAX_VALUE;
164 for (int i = -12; i < -2; ++i) {
165 TestProblem1 pb = new TestProblem1();
166 double minStep = 0;
167 double maxStep = pb.getFinalTime() - pb.getInitialTime();
168 double scalAbsoluteTolerance = Math.pow(10.0, i);
169 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
170
171 EmbeddedRungeKuttaIntegrator integ =
172 new DormandPrince54Integrator(minStep, maxStep,
173 scalAbsoluteTolerance, scalRelativeTolerance);
174 TestProblemHandler handler = new TestProblemHandler(pb, integ);
175 integ.setSafety(0.8);
176 integ.setMaxGrowth(5.0);
177 integ.setMinReduction(0.3);
178 integ.addStepHandler(handler);
179 integ.integrate(pb,
180 pb.getInitialTime(), pb.getInitialState(),
181 pb.getFinalTime(), new double[pb.getDimension()]);
182 assertEquals(0.8, integ.getSafety(), 1.0e-12);
183 assertEquals(5.0, integ.getMaxGrowth(), 1.0e-12);
184 assertEquals(0.3, integ.getMinReduction(), 1.0e-12);
185
186
187
188
189 assertTrue(handler.getMaximalValueError() < (0.7 * scalAbsoluteTolerance));
190 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
191
192 int calls = pb.getCalls();
193 assertTrue(calls <= previousCalls);
194 previousCalls = calls;
195
196 }
197
198 }
199
200 public void testEvents()
201 throws DerivativeException, IntegratorException {
202
203 TestProblem4 pb = new TestProblem4();
204 double minStep = 0;
205 double maxStep = pb.getFinalTime() - pb.getInitialTime();
206 double scalAbsoluteTolerance = 1.0e-8;
207 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
208
209 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
210 scalAbsoluteTolerance,
211 scalRelativeTolerance);
212 TestProblemHandler handler = new TestProblemHandler(pb, integ);
213 integ.addStepHandler(handler);
214 EventHandler[] functions = pb.getEventsHandlers();
215 for (int l = 0; l < functions.length; ++l) {
216 integ.addEventHandler(functions[l],
217 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
218 }
219 assertEquals(functions.length, integ.getEventHandlers().size());
220 integ.integrate(pb,
221 pb.getInitialTime(), pb.getInitialState(),
222 pb.getFinalTime(), new double[pb.getDimension()]);
223
224 assertTrue(handler.getMaximalValueError() < 5.0e-6);
225 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
226 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
227 integ.clearEventHandlers();
228 assertEquals(0, integ.getEventHandlers().size());
229
230 }
231
232 public void testKepler()
233 throws DerivativeException, IntegratorException {
234
235 final TestProblem3 pb = new TestProblem3(0.9);
236 double minStep = 0;
237 double maxStep = pb.getFinalTime() - pb.getInitialTime();
238 double scalAbsoluteTolerance = 1.0e-8;
239 double scalRelativeTolerance = scalAbsoluteTolerance;
240
241 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
242 scalAbsoluteTolerance,
243 scalRelativeTolerance);
244 integ.addStepHandler(new KeplerHandler(pb));
245 integ.integrate(pb,
246 pb.getInitialTime(), pb.getInitialState(),
247 pb.getFinalTime(), new double[pb.getDimension()]);
248
249 assertTrue(pb.getCalls() < 2800);
250
251 }
252
253 public void testVariableSteps()
254 throws DerivativeException, IntegratorException {
255
256 final TestProblem3 pb = new TestProblem3(0.9);
257 double minStep = 0;
258 double maxStep = pb.getFinalTime() - pb.getInitialTime();
259 double scalAbsoluteTolerance = 1.0e-8;
260 double scalRelativeTolerance = scalAbsoluteTolerance;
261
262 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
263 scalAbsoluteTolerance,
264 scalRelativeTolerance);
265 integ.addStepHandler(new VariableHandler());
266 double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
267 pb.getFinalTime(), new double[pb.getDimension()]);
268 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
269 }
270
271 private static class KeplerHandler implements StepHandler {
272 private static final long serialVersionUID = -1645853847806655456L;
273
274 public KeplerHandler(TestProblem3 pb) {
275 this.pb = pb;
276 reset();
277 }
278 public boolean requiresDenseOutput() {
279 return true;
280 }
281 public void reset() {
282 nbSteps = 0;
283 maxError = 0;
284 }
285 public void handleStep(StepInterpolator interpolator,
286 boolean isLast)
287 throws DerivativeException {
288
289 ++nbSteps;
290 for (int a = 1; a < 10; ++a) {
291
292 double prev = interpolator.getPreviousTime();
293 double curr = interpolator.getCurrentTime();
294 double interp = ((10 - a) * prev + a * curr) / 10;
295 interpolator.setInterpolatedTime(interp);
296
297 double[] interpolatedY = interpolator.getInterpolatedState ();
298 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
299 double dx = interpolatedY[0] - theoreticalY[0];
300 double dy = interpolatedY[1] - theoreticalY[1];
301 double error = dx * dx + dy * dy;
302 if (error > maxError) {
303 maxError = error;
304 }
305 }
306 if (isLast) {
307 assertTrue(maxError < 7.0e-10);
308 assertTrue(nbSteps < 400);
309 }
310 }
311 private int nbSteps;
312 private double maxError;
313 private TestProblem3 pb;
314 }
315
316 private static class VariableHandler implements StepHandler {
317 private static final long serialVersionUID = -5196650833828379228L;
318 public VariableHandler() {
319 firstTime = true;
320 minStep = 0;
321 maxStep = 0;
322 }
323 public boolean requiresDenseOutput() {
324 return false;
325 }
326 public void reset() {
327 firstTime = true;
328 minStep = 0;
329 maxStep = 0;
330 }
331 public void handleStep(StepInterpolator interpolator,
332 boolean isLast) {
333
334 double step = Math.abs(interpolator.getCurrentTime()
335 - interpolator.getPreviousTime());
336 if (firstTime) {
337 minStep = Math.abs(step);
338 maxStep = minStep;
339 firstTime = false;
340 } else {
341 if (step < minStep) {
342 minStep = step;
343 }
344 if (step > maxStep) {
345 maxStep = step;
346 }
347 }
348
349 if (isLast) {
350 assertTrue(minStep < (1.0 / 450.0));
351 assertTrue(maxStep > (1.0 / 4.2));
352 }
353 }
354 private boolean firstTime;
355 private double minStep;
356 private double maxStep;
357 }
358
359 public static Test suite() {
360 return new TestSuite(DormandPrince54IntegratorTest.class);
361 }
362
363 }