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