001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode.nonstiff;
019    
020    import org.apache.commons.math.ode.DerivativeException;
021    import org.apache.commons.math.ode.FirstOrderIntegrator;
022    import org.apache.commons.math.ode.IntegratorException;
023    import org.apache.commons.math.ode.TestProblem1;
024    import org.apache.commons.math.ode.TestProblem3;
025    import org.apache.commons.math.ode.TestProblem4;
026    import org.apache.commons.math.ode.TestProblem5;
027    import org.apache.commons.math.ode.TestProblemAbstract;
028    import org.apache.commons.math.ode.TestProblemHandler;
029    import org.apache.commons.math.ode.events.EventHandler;
030    import org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator;
031    import org.apache.commons.math.ode.sampling.StepHandler;
032    import org.apache.commons.math.ode.sampling.StepInterpolator;
033    
034    import junit.framework.*;
035    
036    public class GraggBulirschStoerIntegratorTest
037      extends TestCase {
038    
039      public GraggBulirschStoerIntegratorTest(String name) {
040        super(name);
041      }
042    
043      public void testDimensionCheck() {
044        try  {
045          TestProblem1 pb = new TestProblem1();
046          AdaptiveStepsizeIntegrator integrator =
047            new GraggBulirschStoerIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
048          integrator.integrate(pb,
049                               0.0, new double[pb.getDimension()+10],
050                               1.0, new double[pb.getDimension()+10]);
051          fail("an exception should have been thrown");
052        } catch(DerivativeException de) {
053          fail("wrong exception caught");
054        } catch(IntegratorException ie) {
055        }
056      }
057    
058      public void testNullIntervalCheck() {
059        try  {
060          TestProblem1 pb = new TestProblem1();
061          GraggBulirschStoerIntegrator integrator =
062            new GraggBulirschStoerIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10);
063          integrator.integrate(pb,
064                               0.0, new double[pb.getDimension()],
065                               0.0, new double[pb.getDimension()]);
066          fail("an exception should have been thrown");
067        } catch(DerivativeException de) {
068          fail("wrong exception caught");
069        } catch(IntegratorException ie) {
070        }
071      }
072    
073      public void testMinStep() {
074    
075        try {
076          TestProblem5 pb  = new TestProblem5();
077          double minStep   = 0.1 * Math.abs(pb.getFinalTime() - pb.getInitialTime());
078          double maxStep   = Math.abs(pb.getFinalTime() - pb.getInitialTime());
079          double[] vecAbsoluteTolerance = { 1.0e-20, 1.0e-21 };
080          double[] vecRelativeTolerance = { 1.0e-20, 1.0e-21 };
081    
082          FirstOrderIntegrator integ =
083            new GraggBulirschStoerIntegrator(minStep, maxStep,
084                                             vecAbsoluteTolerance, vecRelativeTolerance);
085          TestProblemHandler handler = new TestProblemHandler(pb, integ);
086          integ.addStepHandler(handler);
087          integ.integrate(pb,
088                          pb.getInitialTime(), pb.getInitialState(),
089                          pb.getFinalTime(), new double[pb.getDimension()]);
090          fail("an exception should have been thrown");
091        } catch(DerivativeException de) {
092          fail("wrong exception caught");
093        } catch(IntegratorException ie) {
094        }
095    
096      }
097    
098      public void testBackward()
099          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    }