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