1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
21 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22 import org.apache.commons.math4.legacy.exception.NoBracketingException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24 import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
25 import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
26 import org.apache.commons.math4.legacy.ode.TestProblem1;
27 import org.apache.commons.math4.legacy.ode.TestProblem3;
28 import org.apache.commons.math4.legacy.ode.TestProblem4;
29 import org.apache.commons.math4.legacy.ode.TestProblem5;
30 import org.apache.commons.math4.legacy.ode.TestProblemAbstract;
31 import org.apache.commons.math4.legacy.ode.TestProblemHandler;
32 import org.apache.commons.math4.legacy.ode.events.EventHandler;
33 import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
34 import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
35 import org.apache.commons.math4.core.jdkmath.JdkMath;
36 import org.junit.Assert;
37 import org.junit.Test;
38
39
40 public class GraggBulirschStoerIntegratorTest {
41
42 @Test(expected=DimensionMismatchException.class)
43 public void testDimensionCheck()
44 throws DimensionMismatchException, NumberIsTooSmallException,
45 MaxCountExceededException, NoBracketingException {
46 TestProblem1 pb = new TestProblem1();
47 AdaptiveStepsizeIntegrator integrator =
48 new GraggBulirschStoerIntegrator(0.0, 1.0, 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 }
53
54 @Test(expected=NumberIsTooSmallException.class)
55 public void testNullIntervalCheck()
56 throws DimensionMismatchException, NumberIsTooSmallException,
57 MaxCountExceededException, NoBracketingException {
58 TestProblem1 pb = new TestProblem1();
59 GraggBulirschStoerIntegrator integrator =
60 new GraggBulirschStoerIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
61 integrator.integrate(pb,
62 0.0, new double[pb.getDimension()],
63 0.0, new double[pb.getDimension()]);
64 }
65
66 @Test(expected=NumberIsTooSmallException.class)
67 public void testMinStep()
68 throws DimensionMismatchException, NumberIsTooSmallException,
69 MaxCountExceededException, NoBracketingException {
70
71 TestProblem5 pb = new TestProblem5();
72 double minStep = 0.1 * JdkMath.abs(pb.getFinalTime() - pb.getInitialTime());
73 double maxStep = JdkMath.abs(pb.getFinalTime() - pb.getInitialTime());
74 double[] vecAbsoluteTolerance = { 1.0e-20, 1.0e-21 };
75 double[] vecRelativeTolerance = { 1.0e-20, 1.0e-21 };
76
77 FirstOrderIntegrator integ =
78 new GraggBulirschStoerIntegrator(minStep, maxStep,
79 vecAbsoluteTolerance, vecRelativeTolerance);
80 TestProblemHandler handler = new TestProblemHandler(pb, integ);
81 integ.addStepHandler(handler);
82 integ.integrate(pb,
83 pb.getInitialTime(), pb.getInitialState(),
84 pb.getFinalTime(), new double[pb.getDimension()]);
85 }
86
87 @Test
88 public void testBackward()
89 throws DimensionMismatchException, NumberIsTooSmallException,
90 MaxCountExceededException, NoBracketingException {
91
92 TestProblem5 pb = new TestProblem5();
93 double minStep = 0;
94 double maxStep = pb.getFinalTime() - pb.getInitialTime();
95 double scalAbsoluteTolerance = 1.0e-8;
96 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
97
98 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep,
99 scalAbsoluteTolerance,
100 scalRelativeTolerance);
101 TestProblemHandler handler = new TestProblemHandler(pb, integ);
102 integ.addStepHandler(handler);
103 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
104 pb.getFinalTime(), new double[pb.getDimension()]);
105
106 Assert.assertTrue(handler.getLastError() < 7.5e-9);
107 Assert.assertTrue(handler.getMaximalValueError() < 8.1e-9);
108 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
109 Assert.assertEquals("Gragg-Bulirsch-Stoer", integ.getName());
110 }
111
112 @Test
113 public void testIncreasingTolerance()
114 throws DimensionMismatchException, NumberIsTooSmallException,
115 MaxCountExceededException, NoBracketingException {
116
117 int previousCalls = Integer.MAX_VALUE;
118 for (int i = -12; i < -4; ++i) {
119 TestProblem1 pb = new TestProblem1();
120 double minStep = 0;
121 double maxStep = pb.getFinalTime() - pb.getInitialTime();
122 double absTolerance = JdkMath.pow(10.0, i);
123 double relTolerance = absTolerance;
124
125 FirstOrderIntegrator integ =
126 new GraggBulirschStoerIntegrator(minStep, maxStep,
127 absTolerance, relTolerance);
128 TestProblemHandler handler = new TestProblemHandler(pb, integ);
129 integ.addStepHandler(handler);
130 integ.integrate(pb,
131 pb.getInitialTime(), pb.getInitialState(),
132 pb.getFinalTime(), new double[pb.getDimension()]);
133
134
135
136
137 double ratio = handler.getMaximalValueError() / absTolerance;
138 Assert.assertTrue(ratio < 2.4);
139 Assert.assertTrue(ratio > 0.02);
140 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
141
142 int calls = pb.getCalls();
143 Assert.assertEquals(integ.getEvaluations(), calls);
144 Assert.assertTrue(calls <= previousCalls);
145 previousCalls = calls;
146 }
147 }
148
149 @Test
150 public void testIntegratorControls()
151 throws DimensionMismatchException, NumberIsTooSmallException,
152 MaxCountExceededException, NoBracketingException {
153
154 TestProblem3 pb = new TestProblem3(0.999);
155 GraggBulirschStoerIntegrator integ =
156 new GraggBulirschStoerIntegrator(0, pb.getFinalTime() - pb.getInitialTime(),
157 1.0e-8, 1.0e-10);
158
159 double errorWithDefaultSettings = getMaxError(integ, pb);
160
161
162 integ.setStabilityCheck(true, 2, 1, 0.99);
163 Assert.assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
164 integ.setStabilityCheck(true, -1, -1, -1);
165
166 integ.setControlFactors(0.5, 0.99, 0.1, 2.5);
167 Assert.assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
168 integ.setControlFactors(-1, -1, -1, -1);
169
170 integ.setOrderControl(10, 0.7, 0.95);
171 Assert.assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
172 integ.setOrderControl(-1, -1, -1);
173
174 integ.setInterpolationControl(true, 3);
175 Assert.assertTrue(errorWithDefaultSettings < getMaxError(integ, pb));
176 integ.setInterpolationControl(true, -1);
177 }
178
179 private double getMaxError(FirstOrderIntegrator integrator, TestProblemAbstract pb)
180 throws DimensionMismatchException, NumberIsTooSmallException,
181 MaxCountExceededException, NoBracketingException {
182 TestProblemHandler handler = new TestProblemHandler(pb, integrator);
183 integrator.addStepHandler(handler);
184 integrator.integrate(pb,
185 pb.getInitialTime(), pb.getInitialState(),
186 pb.getFinalTime(), new double[pb.getDimension()]);
187 return handler.getMaximalValueError();
188 }
189
190 @Test
191 public void testEvents()
192 throws DimensionMismatchException, NumberIsTooSmallException,
193 MaxCountExceededException, NoBracketingException {
194
195 TestProblem4 pb = new TestProblem4();
196 double minStep = 0;
197 double maxStep = pb.getFinalTime() - pb.getInitialTime();
198 double scalAbsoluteTolerance = 1.0e-10;
199 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
200
201 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep,
202 scalAbsoluteTolerance,
203 scalRelativeTolerance);
204 TestProblemHandler handler = new TestProblemHandler(pb, integ);
205 integ.addStepHandler(handler);
206 EventHandler[] functions = pb.getEventsHandlers();
207 double convergence = 1.0e-8 * maxStep;
208 for (int l = 0; l < functions.length; ++l) {
209 integ.addEventHandler(functions[l], Double.POSITIVE_INFINITY, convergence, 1000);
210 }
211 Assert.assertEquals(functions.length, integ.getEventHandlers().size());
212 integ.integrate(pb,
213 pb.getInitialTime(), pb.getInitialState(),
214 pb.getFinalTime(), new double[pb.getDimension()]);
215
216 Assert.assertTrue(handler.getMaximalValueError() < 4.0e-7);
217 Assert.assertEquals(0, handler.getMaximalTimeError(), convergence);
218 Assert.assertEquals(12.0, handler.getLastTime(), convergence);
219 integ.clearEventHandlers();
220 Assert.assertEquals(0, integ.getEventHandlers().size());
221 }
222
223 @Test
224 public void testKepler()
225 throws DimensionMismatchException, NumberIsTooSmallException,
226 MaxCountExceededException, NoBracketingException {
227
228 final TestProblem3 pb = new TestProblem3(0.9);
229 double minStep = 0;
230 double maxStep = pb.getFinalTime() - pb.getInitialTime();
231 double absTolerance = 1.0e-6;
232 double relTolerance = 1.0e-6;
233
234 FirstOrderIntegrator integ =
235 new GraggBulirschStoerIntegrator(minStep, maxStep,
236 absTolerance, relTolerance);
237 integ.addStepHandler(new KeplerStepHandler(pb));
238 integ.integrate(pb,
239 pb.getInitialTime(), pb.getInitialState(),
240 pb.getFinalTime(), new double[pb.getDimension()]);
241
242 Assert.assertEquals(integ.getEvaluations(), pb.getCalls());
243 Assert.assertTrue(pb.getCalls() < 2150);
244 }
245
246 @Test
247 public void testVariableSteps()
248 throws DimensionMismatchException, NumberIsTooSmallException,
249 MaxCountExceededException, NoBracketingException {
250
251 final TestProblem3 pb = new TestProblem3(0.9);
252 double minStep = 0;
253 double maxStep = pb.getFinalTime() - pb.getInitialTime();
254 double absTolerance = 1.0e-8;
255 double relTolerance = 1.0e-8;
256 FirstOrderIntegrator integ =
257 new GraggBulirschStoerIntegrator(minStep, maxStep,
258 absTolerance, relTolerance);
259 integ.addStepHandler(new VariableStepHandler());
260 double stopTime = integ.integrate(pb,
261 pb.getInitialTime(), pb.getInitialState(),
262 pb.getFinalTime(), new double[pb.getDimension()]);
263 Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
264 Assert.assertEquals("Gragg-Bulirsch-Stoer", integ.getName());
265 }
266
267 @Test
268 public void testTooLargeFirstStep()
269 throws DimensionMismatchException, NumberIsTooSmallException,
270 MaxCountExceededException, NoBracketingException {
271
272 AdaptiveStepsizeIntegrator integ =
273 new GraggBulirschStoerIntegrator(0, Double.POSITIVE_INFINITY, Double.NaN, Double.NaN);
274 final double start = 0.0;
275 final double end = 0.001;
276 FirstOrderDifferentialEquations equations = new FirstOrderDifferentialEquations() {
277
278 @Override
279 public int getDimension() {
280 return 1;
281 }
282
283 @Override
284 public void computeDerivatives(double t, double[] y, double[] yDot) {
285 Assert.assertTrue(t >= JdkMath.nextAfter(start, Double.NEGATIVE_INFINITY));
286 Assert.assertTrue(t <= JdkMath.nextAfter(end, Double.POSITIVE_INFINITY));
287 yDot[0] = -100.0 * y[0];
288 }
289 };
290
291 integ.setStepSizeControl(0, 1.0, 1.0e-6, 1.0e-8);
292 integ.integrate(equations, start, new double[] { 1.0 }, end, new double[1]);
293 }
294
295 @Test
296 public void testUnstableDerivative()
297 throws DimensionMismatchException, NumberIsTooSmallException,
298 MaxCountExceededException, NoBracketingException {
299 final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
300 FirstOrderIntegrator integ =
301 new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
302 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
303 double[] y = { Double.NaN };
304 integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
305 Assert.assertEquals(8.0, y[0], 1.0e-12);
306 }
307
308 @Test
309 public void testIssue596()
310 throws DimensionMismatchException, NumberIsTooSmallException,
311 MaxCountExceededException, NoBracketingException {
312 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(1e-10, 100.0, 1e-7, 1e-7);
313 integ.addStepHandler(new StepHandler() {
314
315 @Override
316 public void init(double t0, double[] y0, double t) {
317 }
318
319 @Override
320 public void handleStep(StepInterpolator interpolator, boolean isLast)
321 throws MaxCountExceededException {
322 double t = interpolator.getCurrentTime();
323 interpolator.setInterpolatedTime(t);
324 double[] y = interpolator.getInterpolatedState();
325 double[] yDot = interpolator.getInterpolatedDerivatives();
326 Assert.assertEquals(3.0 * t - 5.0, y[0], 1.0e-14);
327 Assert.assertEquals(3.0, yDot[0], 1.0e-14);
328 }
329 });
330 double[] y = {4.0};
331 double t0 = 3.0;
332 double tend = 10.0;
333 integ.integrate(new FirstOrderDifferentialEquations() {
334 @Override
335 public int getDimension() {
336 return 1;
337 }
338
339 @Override
340 public void computeDerivatives(double t, double[] y, double[] yDot) {
341 yDot[0] = 3.0;
342 }
343 }, t0, y, tend, y);
344 }
345
346 private static final class KeplerStepHandler implements StepHandler {
347 KeplerStepHandler(TestProblem3 pb) {
348 this.pb = pb;
349 }
350 @Override
351 public void init(double t0, double[] y0, double t) {
352 nbSteps = 0;
353 maxError = 0;
354 }
355 @Override
356 public void handleStep(StepInterpolator interpolator, boolean isLast)
357 throws MaxCountExceededException {
358
359 ++nbSteps;
360 for (int a = 1; a < 100; ++a) {
361
362 double prev = interpolator.getPreviousTime();
363 double curr = interpolator.getCurrentTime();
364 double interp = ((100 - a) * prev + a * curr) / 100;
365 interpolator.setInterpolatedTime(interp);
366
367 double[] interpolatedY = interpolator.getInterpolatedState ();
368 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
369 double dx = interpolatedY[0] - theoreticalY[0];
370 double dy = interpolatedY[1] - theoreticalY[1];
371 double error = dx * dx + dy * dy;
372 if (error > maxError) {
373 maxError = error;
374 }
375 }
376 if (isLast) {
377 Assert.assertTrue(maxError < 2.7e-6);
378 Assert.assertTrue(nbSteps < 80);
379 }
380 }
381 private int nbSteps;
382 private double maxError;
383 private TestProblem3 pb;
384 }
385
386 public static class VariableStepHandler implements StepHandler {
387 public VariableStepHandler() {
388 firstTime = true;
389 minStep = 0;
390 maxStep = 0;
391 }
392 @Override
393 public void init(double t0, double[] y0, double t) {
394 firstTime = true;
395 minStep = 0;
396 maxStep = 0;
397 }
398 @Override
399 public void handleStep(StepInterpolator interpolator,
400 boolean isLast) {
401
402 double step = JdkMath.abs(interpolator.getCurrentTime()
403 - interpolator.getPreviousTime());
404 if (firstTime) {
405 minStep = JdkMath.abs(step);
406 maxStep = minStep;
407 firstTime = false;
408 } else {
409 if (step < minStep) {
410 minStep = step;
411 }
412 if (step > maxStep) {
413 maxStep = step;
414 }
415 }
416
417 if (isLast) {
418 Assert.assertTrue(minStep < 8.2e-3);
419 Assert.assertTrue(maxStep > 1.5);
420 }
421 }
422 private boolean firstTime;
423 private double minStep;
424 private double maxStep;
425 }
426 }