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