1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode;
19
20 import org.apache.commons.math.ode.AdaptiveStepsizeIntegrator;
21 import org.apache.commons.math.ode.DerivativeException;
22 import org.apache.commons.math.ode.DormandPrince54Integrator;
23 import org.apache.commons.math.ode.FirstOrderIntegrator;
24 import org.apache.commons.math.ode.IntegratorException;
25 import org.apache.commons.math.ode.StepHandler;
26 import org.apache.commons.math.ode.StepInterpolator;
27 import org.apache.commons.math.ode.SwitchingFunction;
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.setStepHandler(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.setStepHandler(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 private static class DP54SmallLastHandler implements StepHandler {
105
106 public DP54SmallLastHandler(double minStep) {
107 lastSeen = false;
108 this.minStep = minStep;
109 }
110
111 public boolean requiresDenseOutput() {
112 return false;
113 }
114
115 public void reset() {
116 }
117
118 public void handleStep(StepInterpolator interpolator, boolean isLast) {
119 if (isLast) {
120 lastSeen = true;
121 double h = interpolator.getCurrentTime() - interpolator.getPreviousTime();
122 assertTrue(Math.abs(h) < minStep);
123 }
124 }
125
126 public boolean wasLastSeen() {
127 return lastSeen;
128 }
129
130 private boolean lastSeen;
131 private double minStep;
132
133 }
134
135 public void testIncreasingTolerance()
136 throws DerivativeException, IntegratorException {
137
138 int previousCalls = Integer.MAX_VALUE;
139 for (int i = -12; i < -2; ++i) {
140 TestProblem1 pb = new TestProblem1();
141 double minStep = 0;
142 double maxStep = pb.getFinalTime() - pb.getInitialTime();
143 double scalAbsoluteTolerance = Math.pow(10.0, i);
144 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
145
146 EmbeddedRungeKuttaIntegrator integ =
147 new DormandPrince54Integrator(minStep, maxStep,
148 scalAbsoluteTolerance, scalRelativeTolerance);
149 TestProblemHandler handler = new TestProblemHandler(pb, integ);
150 integ.setSafety(0.8);
151 integ.setMaxGrowth(5.0);
152 integ.setMinReduction(0.3);
153 integ.setStepHandler(handler);
154 integ.integrate(pb,
155 pb.getInitialTime(), pb.getInitialState(),
156 pb.getFinalTime(), new double[pb.getDimension()]);
157 assertEquals(0.8, integ.getSafety(), 1.0e-12);
158 assertEquals(5.0, integ.getMaxGrowth(), 1.0e-12);
159 assertEquals(0.3, integ.getMinReduction(), 1.0e-12);
160
161
162
163
164 assertTrue(handler.getMaximalValueError() < (0.7 * scalAbsoluteTolerance));
165 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
166
167 int calls = pb.getCalls();
168 assertTrue(calls <= previousCalls);
169 previousCalls = calls;
170
171 }
172
173 }
174
175 public void testSwitchingFunctions()
176 throws DerivativeException, IntegratorException {
177
178 TestProblem4 pb = new TestProblem4();
179 double minStep = 0;
180 double maxStep = pb.getFinalTime() - pb.getInitialTime();
181 double scalAbsoluteTolerance = 1.0e-8;
182 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
183
184 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
185 scalAbsoluteTolerance,
186 scalRelativeTolerance);
187 TestProblemHandler handler = new TestProblemHandler(pb, integ);
188 integ.setStepHandler(handler);
189 SwitchingFunction[] functions = pb.getSwitchingFunctions();
190 for (int l = 0; l < functions.length; ++l) {
191 integ.addSwitchingFunction(functions[l],
192 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
193 }
194 integ.integrate(pb,
195 pb.getInitialTime(), pb.getInitialState(),
196 pb.getFinalTime(), new double[pb.getDimension()]);
197
198 assertTrue(handler.getMaximalValueError() < 5.0e-6);
199 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
200 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
201
202 }
203
204 public void testKepler()
205 throws DerivativeException, IntegratorException {
206
207 final TestProblem3 pb = new TestProblem3(0.9);
208 double minStep = 0;
209 double maxStep = pb.getFinalTime() - pb.getInitialTime();
210 double scalAbsoluteTolerance = 1.0e-8;
211 double scalRelativeTolerance = scalAbsoluteTolerance;
212
213 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
214 scalAbsoluteTolerance,
215 scalRelativeTolerance);
216 integ.setStepHandler(new KeplerHandler(pb));
217 integ.integrate(pb,
218 pb.getInitialTime(), pb.getInitialState(),
219 pb.getFinalTime(), new double[pb.getDimension()]);
220
221 assertTrue(pb.getCalls() < 2800);
222
223 }
224
225 public void testVariableSteps()
226 throws DerivativeException, IntegratorException {
227
228 final TestProblem3 pb = new TestProblem3(0.9);
229 double minStep = 0;
230 double maxStep = pb.getFinalTime() - pb.getInitialTime();
231 double scalAbsoluteTolerance = 1.0e-8;
232 double scalRelativeTolerance = scalAbsoluteTolerance;
233
234 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
235 scalAbsoluteTolerance,
236 scalRelativeTolerance);
237 integ.setStepHandler(new VariableHandler());
238 integ.integrate(pb,
239 pb.getInitialTime(), pb.getInitialState(),
240 pb.getFinalTime(), new double[pb.getDimension()]);
241 }
242
243 private static class KeplerHandler implements StepHandler {
244 public KeplerHandler(TestProblem3 pb) {
245 this.pb = pb;
246 reset();
247 }
248 public boolean requiresDenseOutput() {
249 return true;
250 }
251 public void reset() {
252 nbSteps = 0;
253 maxError = 0;
254 }
255 public void handleStep(StepInterpolator interpolator,
256 boolean isLast)
257 throws DerivativeException {
258
259 ++nbSteps;
260 for (int a = 1; a < 10; ++a) {
261
262 double prev = interpolator.getPreviousTime();
263 double curr = interpolator.getCurrentTime();
264 double interp = ((10 - a) * prev + a * curr) / 10;
265 interpolator.setInterpolatedTime(interp);
266
267 double[] interpolatedY = interpolator.getInterpolatedState ();
268 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
269 double dx = interpolatedY[0] - theoreticalY[0];
270 double dy = interpolatedY[1] - theoreticalY[1];
271 double error = dx * dx + dy * dy;
272 if (error > maxError) {
273 maxError = error;
274 }
275 }
276 if (isLast) {
277 assertTrue(maxError < 7.0e-10);
278 assertTrue(nbSteps < 400);
279 }
280 }
281 private int nbSteps;
282 private double maxError;
283 private TestProblem3 pb;
284 }
285
286 private static class VariableHandler implements StepHandler {
287 public VariableHandler() {
288 firstTime = true;
289 minStep = 0;
290 maxStep = 0;
291 }
292 public boolean requiresDenseOutput() {
293 return false;
294 }
295 public void reset() {
296 firstTime = true;
297 minStep = 0;
298 maxStep = 0;
299 }
300 public void handleStep(StepInterpolator interpolator,
301 boolean isLast) {
302
303 double step = Math.abs(interpolator.getCurrentTime()
304 - interpolator.getPreviousTime());
305 if (firstTime) {
306 minStep = Math.abs(step);
307 maxStep = minStep;
308 firstTime = false;
309 } else {
310 if (step < minStep) {
311 minStep = step;
312 }
313 if (step > maxStep) {
314 maxStep = step;
315 }
316 }
317
318 if (isLast) {
319 assertTrue(minStep < (1.0 / 450.0));
320 assertTrue(maxStep > (1.0 / 4.2));
321 }
322 }
323 private boolean firstTime;
324 private double minStep;
325 private double maxStep;
326 }
327
328 public static Test suite() {
329 return new TestSuite(DormandPrince54IntegratorTest.class);
330 }
331
332 }