1   /*
2    * Licensed to the Apache Software Foundation (ASF) 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 ASF 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    *      http://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.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       // the coefficients are only valid for this test
142       // and have been obtained from trial and error
143       // there is no general relation between local and global errors
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     // stability control
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 }