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.GraggBulirschStoerIntegrator;
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 GraggBulirschStoerIntegratorTest
37 extends TestCase {
38
39 public GraggBulirschStoerIntegratorTest(String name) {
40 super(name);
41 }
42
43 public void testDimensionCheck() {
44 try {
45 TestProblem1 pb = new TestProblem1();
46 AdaptiveStepsizeIntegrator integrator =
47 new GraggBulirschStoerIntegrator(0.0, 1.0, 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 GraggBulirschStoerIntegrator integrator =
62 new GraggBulirschStoerIntegrator(0.0, 1.0, 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 TestProblem5 pb = new TestProblem5();
77 double minStep = 0.1 * Math.abs(pb.getFinalTime() - pb.getInitialTime());
78 double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
79 double[] vecAbsoluteTolerance = { 1.0e-20, 1.0e-21 };
80 double[] vecRelativeTolerance = { 1.0e-20, 1.0e-21 };
81
82 FirstOrderIntegrator integ =
83 new GraggBulirschStoerIntegrator(minStep, maxStep,
84 vecAbsoluteTolerance, 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 testBackward()
99 throws DerivativeException, IntegratorException {
100
101 TestProblem5 pb = new TestProblem5();
102 double minStep = 0;
103 double maxStep = pb.getFinalTime() - pb.getInitialTime();
104 double scalAbsoluteTolerance = 1.0e-8;
105 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
106
107 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep,
108 scalAbsoluteTolerance,
109 scalRelativeTolerance);
110 TestProblemHandler handler = new TestProblemHandler(pb, integ);
111 integ.addStepHandler(handler);
112 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
113 pb.getFinalTime(), new double[pb.getDimension()]);
114
115 assertTrue(handler.getLastError() < 9.0e-10);
116 assertTrue(handler.getMaximalValueError() < 9.0e-10);
117 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
118 assertEquals("Gragg-Bulirsch-Stoer", integ.getName());
119 }
120
121 public void testIncreasingTolerance()
122 throws DerivativeException, IntegratorException {
123
124 int previousCalls = Integer.MAX_VALUE;
125 for (int i = -12; i < -4; ++i) {
126 TestProblem1 pb = new TestProblem1();
127 double minStep = 0;
128 double maxStep = pb.getFinalTime() - pb.getInitialTime();
129 double absTolerance = Math.pow(10.0, i);
130 double relTolerance = absTolerance;
131
132 FirstOrderIntegrator integ =
133 new GraggBulirschStoerIntegrator(minStep, maxStep,
134 absTolerance, relTolerance);
135 TestProblemHandler handler = new TestProblemHandler(pb, integ);
136 integ.addStepHandler(handler);
137 integ.integrate(pb,
138 pb.getInitialTime(), pb.getInitialState(),
139 pb.getFinalTime(), new double[pb.getDimension()]);
140
141
142
143
144 double ratio = handler.getMaximalValueError() / absTolerance;
145 assertTrue(ratio < 2.4);
146 assertTrue(ratio > 0.02);
147 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
148
149 int calls = pb.getCalls();
150 assertEquals(integ.getEvaluations(), calls);
151 assertTrue(calls <= previousCalls);
152 previousCalls = calls;
153
154 }
155
156 }
157
158 public void testIntegratorControls()
159 throws DerivativeException, IntegratorException {
160
161 TestProblem3 pb = new TestProblem3(0.999);
162 GraggBulirschStoerIntegrator integ =
163 new GraggBulirschStoerIntegrator(0, pb.getFinalTime() - pb.getInitialTime(),
164 1.0e-8, 1.0e-10);
165
166 double errorWithDefaultSettings = getMaxError(integ, pb);
167
168
169 integ.setStabilityCheck(true, 2, 1, 0.99);
170 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
171 integ.setStabilityCheck(true, -1, -1, -1);
172
173 integ.setStepsizeControl(0.5, 0.99, 0.1, 2.5);
174 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
175 integ.setStepsizeControl(-1, -1, -1, -1);
176
177 integ.setOrderControl(10, 0.7, 0.95);
178 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
179 integ.setOrderControl(-1, -1, -1);
180
181 integ.setInterpolationControl(true, 3);
182 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
183 integ.setInterpolationControl(true, -1);
184
185 }
186
187 private double getMaxError(FirstOrderIntegrator integrator, TestProblemAbstract pb)
188 throws DerivativeException, IntegratorException {
189 TestProblemHandler handler = new TestProblemHandler(pb, integrator);
190 integrator.addStepHandler(handler);
191 integrator.integrate(pb,
192 pb.getInitialTime(), pb.getInitialState(),
193 pb.getFinalTime(), new double[pb.getDimension()]);
194 return handler.getMaximalValueError();
195 }
196
197 public void testEvents()
198 throws DerivativeException, IntegratorException {
199
200 TestProblem4 pb = new TestProblem4();
201 double minStep = 0;
202 double maxStep = pb.getFinalTime() - pb.getInitialTime();
203 double scalAbsoluteTolerance = 1.0e-10;
204 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
205
206 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep,
207 scalAbsoluteTolerance,
208 scalRelativeTolerance);
209 TestProblemHandler handler = new TestProblemHandler(pb, integ);
210 integ.addStepHandler(handler);
211 EventHandler[] functions = pb.getEventsHandlers();
212 for (int l = 0; l < functions.length; ++l) {
213 integ.addEventHandler(functions[l],
214 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
215 }
216 assertEquals(functions.length, integ.getEventHandlers().size());
217 integ.integrate(pb,
218 pb.getInitialTime(), pb.getInitialState(),
219 pb.getFinalTime(), new double[pb.getDimension()]);
220
221 assertTrue(handler.getMaximalValueError() < 5.0e-8);
222 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
223 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
224 integ.clearEventHandlers();
225 assertEquals(0, integ.getEventHandlers().size());
226
227 }
228
229 public void testKepler()
230 throws DerivativeException, IntegratorException {
231
232 final TestProblem3 pb = new TestProblem3(0.9);
233 double minStep = 0;
234 double maxStep = pb.getFinalTime() - pb.getInitialTime();
235 double absTolerance = 1.0e-6;
236 double relTolerance = 1.0e-6;
237
238 FirstOrderIntegrator integ =
239 new GraggBulirschStoerIntegrator(minStep, maxStep,
240 absTolerance, relTolerance);
241 integ.addStepHandler(new KeplerStepHandler(pb));
242 integ.integrate(pb,
243 pb.getInitialTime(), pb.getInitialState(),
244 pb.getFinalTime(), new double[pb.getDimension()]);
245
246 assertEquals(integ.getEvaluations(), pb.getCalls());
247 assertTrue(pb.getCalls() < 2150);
248
249 }
250
251 public void testVariableSteps()
252 throws DerivativeException, IntegratorException {
253
254 final TestProblem3 pb = new TestProblem3(0.9);
255 double minStep = 0;
256 double maxStep = pb.getFinalTime() - pb.getInitialTime();
257 double absTolerance = 1.0e-8;
258 double relTolerance = 1.0e-8;
259 FirstOrderIntegrator integ =
260 new GraggBulirschStoerIntegrator(minStep, maxStep,
261 absTolerance, relTolerance);
262 integ.addStepHandler(new VariableStepHandler());
263 double stopTime = integ.integrate(pb,
264 pb.getInitialTime(), pb.getInitialState(),
265 pb.getFinalTime(), new double[pb.getDimension()]);
266 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
267 assertEquals("Gragg-Bulirsch-Stoer", integ.getName());
268 }
269
270 public void testUnstableDerivative()
271 throws DerivativeException, IntegratorException {
272 final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
273 FirstOrderIntegrator integ =
274 new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
275 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
276 double[] y = { Double.NaN };
277 integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
278 assertEquals(8.0, y[0], 1.0e-12);
279 }
280
281 private static class KeplerStepHandler implements StepHandler {
282 public KeplerStepHandler(TestProblem3 pb) {
283 this.pb = pb;
284 reset();
285 }
286 public boolean requiresDenseOutput() {
287 return true;
288 }
289 public void reset() {
290 nbSteps = 0;
291 maxError = 0;
292 }
293 public void handleStep(StepInterpolator interpolator,
294 boolean isLast)
295 throws DerivativeException {
296
297 ++nbSteps;
298 for (int a = 1; a < 100; ++a) {
299
300 double prev = interpolator.getPreviousTime();
301 double curr = interpolator.getCurrentTime();
302 double interp = ((100 - a) * prev + a * curr) / 100;
303 interpolator.setInterpolatedTime(interp);
304
305 double[] interpolatedY = interpolator.getInterpolatedState ();
306 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
307 double dx = interpolatedY[0] - theoreticalY[0];
308 double dy = interpolatedY[1] - theoreticalY[1];
309 double error = dx * dx + dy * dy;
310 if (error > maxError) {
311 maxError = error;
312 }
313 }
314 if (isLast) {
315 assertTrue(maxError < 2.7e-6);
316 assertTrue(nbSteps < 80);
317 }
318 }
319 private int nbSteps;
320 private double maxError;
321 private TestProblem3 pb;
322 }
323
324 public static class VariableStepHandler implements StepHandler {
325 public VariableStepHandler() {
326 reset();
327 }
328 public boolean requiresDenseOutput() {
329 return false;
330 }
331 public void reset() {
332 firstTime = true;
333 minStep = 0;
334 maxStep = 0;
335 }
336 public void handleStep(StepInterpolator interpolator,
337 boolean isLast) {
338
339 double step = Math.abs(interpolator.getCurrentTime()
340 - interpolator.getPreviousTime());
341 if (firstTime) {
342 minStep = Math.abs(step);
343 maxStep = minStep;
344 firstTime = false;
345 } else {
346 if (step < minStep) {
347 minStep = step;
348 }
349 if (step > maxStep) {
350 maxStep = step;
351 }
352 }
353
354 if (isLast) {
355 assertTrue(minStep < 8.2e-3);
356 assertTrue(maxStep > 1.7);
357 }
358 }
359 private boolean firstTime;
360 private double minStep;
361 private double maxStep;
362 }
363 public static Test suite() {
364 return new TestSuite(GraggBulirschStoerIntegratorTest.class);
365 }
366
367 }