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