View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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             // the coefficients are only valid for this test
125             // and have been obtained from trial and error
126             // there is no general relation between local and global errors
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         // stability control
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         // since state is approx. linear at g=0 need convergence <= (state tolerance) / 2.
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         // integration error builds up by the last event,
202         // so tolerance is slightly more than the convergence.
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 }