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.TestProblemHandler;
028    import org.apache.commons.math.ode.events.EventHandler;
029    import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
030    import org.apache.commons.math.ode.sampling.DummyStepHandler;
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 DormandPrince853IntegratorTest
037      extends TestCase {
038    
039      public DormandPrince853IntegratorTest(String name) {
040        super(name);
041      }
042    
043      public void testDimensionCheck() {
044        try  {
045          TestProblem1 pb = new TestProblem1();
046          DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.0, 1.0,
047                                                                                 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          DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.0, 1.0,
062                                                                                 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          TestProblem1 pb = new TestProblem1();
077          double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
078          double maxStep = pb.getFinalTime() - pb.getInitialTime();
079          double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
080          double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
081    
082          FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
083                                                                      vecAbsoluteTolerance,
084                                                                      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 testIncreasingTolerance()
099        throws DerivativeException, IntegratorException {
100    
101        int previousCalls = Integer.MAX_VALUE;
102        for (int i = -12; i < -2; ++i) {
103          TestProblem1 pb = new TestProblem1();
104          double minStep = 0;
105          double maxStep = pb.getFinalTime() - pb.getInitialTime();
106          double scalAbsoluteTolerance = Math.pow(10.0, i);
107          double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
108    
109          FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
110                                                                      scalAbsoluteTolerance,
111                                                                      scalRelativeTolerance);
112          TestProblemHandler handler = new TestProblemHandler(pb, integ);
113          integ.addStepHandler(handler);
114          integ.integrate(pb,
115                          pb.getInitialTime(), pb.getInitialState(),
116                          pb.getFinalTime(), new double[pb.getDimension()]);
117    
118          // the 1.3 factor is only valid for this test
119          // and has been obtained from trial and error
120          // there is no general relation between local and global errors
121          assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
122          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
123    
124          int calls = pb.getCalls();
125          assertEquals(integ.getEvaluations(), calls);
126          assertTrue(calls <= previousCalls);
127          previousCalls = calls;
128    
129        }
130    
131      }
132    
133      public void testBackward()
134          throws DerivativeException, IntegratorException {
135    
136          TestProblem5 pb = new TestProblem5();
137          double minStep = 0;
138          double maxStep = pb.getFinalTime() - pb.getInitialTime();
139          double scalAbsoluteTolerance = 1.0e-8;
140          double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
141    
142          FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
143                                                                      scalAbsoluteTolerance,
144                                                                      scalRelativeTolerance);
145          TestProblemHandler handler = new TestProblemHandler(pb, integ);
146          integ.addStepHandler(handler);
147          integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
148                          pb.getFinalTime(), new double[pb.getDimension()]);
149    
150          assertTrue(handler.getLastError() < 8.0e-8);
151          assertTrue(handler.getMaximalValueError() < 2.0e-7);
152          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
153          assertEquals("Dormand-Prince 8 (5, 3)", integ.getName());
154      }
155    
156      public void testEvents()
157        throws DerivativeException, IntegratorException {
158    
159        TestProblem4 pb = new TestProblem4();
160        double minStep = 0;
161        double maxStep = pb.getFinalTime() - pb.getInitialTime();
162        double scalAbsoluteTolerance = 1.0e-9;
163        double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
164    
165        FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
166                                                                    scalAbsoluteTolerance,
167                                                                    scalRelativeTolerance);
168        TestProblemHandler handler = new TestProblemHandler(pb, integ);
169        integ.addStepHandler(handler);
170        EventHandler[] functions = pb.getEventsHandlers();
171        for (int l = 0; l < functions.length; ++l) {
172          integ.addEventHandler(functions[l],
173                                     Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
174        }
175        assertEquals(functions.length, integ.getEventHandlers().size());
176        integ.integrate(pb,
177                        pb.getInitialTime(), pb.getInitialState(),
178                        pb.getFinalTime(), new double[pb.getDimension()]);
179    
180        assertTrue(handler.getMaximalValueError() < 5.0e-8);
181        assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
182        assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
183        integ.clearEventHandlers();
184        assertEquals(0, integ.getEventHandlers().size());
185    
186      }
187    
188      public void testKepler()
189        throws DerivativeException, IntegratorException {
190    
191        final TestProblem3 pb  = new TestProblem3(0.9);
192        double minStep = 0;
193        double maxStep = pb.getFinalTime() - pb.getInitialTime();
194        double scalAbsoluteTolerance = 1.0e-8;
195        double scalRelativeTolerance = scalAbsoluteTolerance;
196    
197        FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
198                                                                    scalAbsoluteTolerance,
199                                                                    scalRelativeTolerance);
200        integ.addStepHandler(new KeplerHandler(pb));
201        integ.integrate(pb,
202                        pb.getInitialTime(), pb.getInitialState(),
203                        pb.getFinalTime(), new double[pb.getDimension()]);
204    
205        assertEquals(integ.getEvaluations(), pb.getCalls());
206        assertTrue(pb.getCalls() < 3300);
207    
208      }
209    
210      public void testVariableSteps()
211        throws DerivativeException, IntegratorException {
212    
213        final TestProblem3 pb  = new TestProblem3(0.9);
214        double minStep = 0;
215        double maxStep = pb.getFinalTime() - pb.getInitialTime();
216        double scalAbsoluteTolerance = 1.0e-8;
217        double scalRelativeTolerance = scalAbsoluteTolerance;
218    
219        FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
220                                                                   scalAbsoluteTolerance,
221                                                                   scalRelativeTolerance);
222        integ.addStepHandler(new VariableHandler());
223        double stopTime = integ.integrate(pb,
224                                          pb.getInitialTime(), pb.getInitialState(),
225                                          pb.getFinalTime(), new double[pb.getDimension()]);
226        assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
227        assertEquals("Dormand-Prince 8 (5, 3)", integ.getName());
228      }
229    
230      public void testNoDenseOutput()
231        throws DerivativeException, IntegratorException {
232        TestProblem1 pb1 = new TestProblem1();
233        TestProblem1 pb2 = pb1.copy();
234        double minStep = 0.1 * (pb1.getFinalTime() - pb1.getInitialTime());
235        double maxStep = pb1.getFinalTime() - pb1.getInitialTime();
236        double scalAbsoluteTolerance = 1.0e-4;
237        double scalRelativeTolerance = 1.0e-4;
238    
239        FirstOrderIntegrator integ = new DormandPrince853Integrator(minStep, maxStep,
240                                                                    scalAbsoluteTolerance,
241                                                                    scalRelativeTolerance);
242        integ.addStepHandler(DummyStepHandler.getInstance());
243        integ.integrate(pb1,
244                        pb1.getInitialTime(), pb1.getInitialState(),
245                        pb1.getFinalTime(), new double[pb1.getDimension()]);
246        int callsWithoutDenseOutput = pb1.getCalls();
247        assertEquals(integ.getEvaluations(), callsWithoutDenseOutput);
248    
249        integ.addStepHandler(new InterpolatingStepHandler());
250        integ.integrate(pb2,
251                        pb2.getInitialTime(), pb2.getInitialState(),
252                        pb2.getFinalTime(), new double[pb2.getDimension()]);
253        int callsWithDenseOutput = pb2.getCalls();
254        assertEquals(integ.getEvaluations(), callsWithDenseOutput);
255    
256        assertTrue(callsWithDenseOutput > callsWithoutDenseOutput);
257    
258      }
259    
260      public void testUnstableDerivative()
261      throws DerivativeException, IntegratorException {
262        final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
263        FirstOrderIntegrator integ =
264          new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
265        integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
266        double[] y = { Double.NaN };
267        integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
268        assertEquals(8.0, y[0], 1.0e-12);
269      }
270    
271      private static class KeplerHandler implements StepHandler {
272        public KeplerHandler(TestProblem3 pb) {
273          this.pb = pb;
274          reset();
275        }
276        public boolean requiresDenseOutput() {
277          return true;
278        }
279        public void reset() {
280          nbSteps = 0;
281          maxError = 0;
282        }
283        public void handleStep(StepInterpolator interpolator,
284                               boolean isLast)
285        throws DerivativeException {
286    
287          ++nbSteps;
288          for (int a = 1; a < 10; ++a) {
289    
290            double prev   = interpolator.getPreviousTime();
291            double curr   = interpolator.getCurrentTime();
292            double interp = ((10 - a) * prev + a * curr) / 10;
293            interpolator.setInterpolatedTime(interp);
294    
295            double[] interpolatedY = interpolator.getInterpolatedState ();
296            double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
297            double dx = interpolatedY[0] - theoreticalY[0];
298            double dy = interpolatedY[1] - theoreticalY[1];
299            double error = dx * dx + dy * dy;
300            if (error > maxError) {
301              maxError = error;
302            }
303          }
304          if (isLast) {
305            assertTrue(maxError < 2.4e-10);
306            assertTrue(nbSteps < 150);
307          }
308        }
309        private int nbSteps;
310        private double maxError;
311        private TestProblem3 pb;
312      }
313    
314      private static class VariableHandler implements StepHandler {
315        public VariableHandler() {
316          reset();
317        }
318        public boolean requiresDenseOutput() {
319          return false;
320        }
321        public void reset() {
322          firstTime = true;
323          minStep = 0;
324          maxStep = 0;
325        }
326        public void handleStep(StepInterpolator interpolator,
327                               boolean isLast) {
328    
329          double step = Math.abs(interpolator.getCurrentTime()
330                                 - interpolator.getPreviousTime());
331          if (firstTime) {
332            minStep   = Math.abs(step);
333            maxStep   = minStep;
334            firstTime = false;
335          } else {
336            if (step < minStep) {
337              minStep = step;
338            }
339            if (step > maxStep) {
340              maxStep = step;
341            }
342          }
343    
344          if (isLast) {
345            assertTrue(minStep < (1.0 / 100.0));
346            assertTrue(maxStep > (1.0 / 2.0));
347          }
348        }
349        private boolean firstTime = true;
350        private double  minStep = 0;
351        private double  maxStep = 0;
352      }
353    
354      private static class InterpolatingStepHandler implements StepHandler {
355        public boolean requiresDenseOutput() {
356          return true;
357        }
358        public void reset() {
359        }
360        public void handleStep(StepInterpolator interpolator,
361                               boolean isLast)
362        throws DerivativeException {
363          double prev = interpolator.getPreviousTime();
364          double curr = interpolator.getCurrentTime();
365          interpolator.setInterpolatedTime(0.5*(prev + curr));
366        }
367      }
368    
369      public static Test suite() {
370        return new TestSuite(DormandPrince853IntegratorTest.class);
371      }
372    
373    }