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.DormandPrince54Integrator;
031    import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
032    import org.apache.commons.math.ode.sampling.StepHandler;
033    import org.apache.commons.math.ode.sampling.StepInterpolator;
034    
035    import junit.framework.*;
036    
037    public class DormandPrince54IntegratorTest
038      extends TestCase {
039    
040      public DormandPrince54IntegratorTest(String name) {
041        super(name);
042      }
043    
044      public void testDimensionCheck() {
045        try  {
046          TestProblem1 pb = new TestProblem1();
047          DormandPrince54Integrator integrator = new DormandPrince54Integrator(0.0, 1.0,
048                                                                               1.0e-10, 1.0e-10);
049          integrator.integrate(pb,
050                               0.0, new double[pb.getDimension()+10],
051                               1.0, new double[pb.getDimension()+10]);
052          fail("an exception should have been thrown");
053        } catch(DerivativeException de) {
054          fail("wrong exception caught");
055        } catch(IntegratorException ie) {
056        }
057      }
058    
059      public void testMinStep() {
060    
061        try {
062          TestProblem1 pb = new TestProblem1();
063          double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
064          double maxStep = pb.getFinalTime() - pb.getInitialTime();
065          double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
066          double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
067    
068          FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
069                                                                     vecAbsoluteTolerance,
070                                                                     vecRelativeTolerance);
071          TestProblemHandler handler = new TestProblemHandler(pb, integ);
072          integ.addStepHandler(handler);
073          integ.integrate(pb,
074                          pb.getInitialTime(), pb.getInitialState(),
075                          pb.getFinalTime(), new double[pb.getDimension()]);
076          fail("an exception should have been thrown");
077        } catch(DerivativeException de) {
078          fail("wrong exception caught");
079        } catch(IntegratorException ie) {
080        }
081    
082      }
083    
084      public void testSmallLastStep()
085        throws DerivativeException, IntegratorException {
086    
087        TestProblemAbstract pb = new TestProblem5();
088        double minStep = 1.25;
089        double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
090        double scalAbsoluteTolerance = 6.0e-4;
091        double scalRelativeTolerance = 6.0e-4;
092    
093        AdaptiveStepsizeIntegrator integ =
094          new DormandPrince54Integrator(minStep, maxStep,
095                                        scalAbsoluteTolerance,
096                                        scalRelativeTolerance);
097    
098        DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
099        integ.addStepHandler(handler);
100        integ.setInitialStepSize(1.7);
101        integ.integrate(pb,
102                        pb.getInitialTime(), pb.getInitialState(),
103                        pb.getFinalTime(), new double[pb.getDimension()]);
104        assertTrue(handler.wasLastSeen());
105        assertEquals("Dormand-Prince 5(4)", integ.getName());
106    
107      }
108    
109      public void testBackward()
110          throws DerivativeException, IntegratorException {
111    
112          TestProblem5 pb = new TestProblem5();
113          double minStep = 0;
114          double maxStep = pb.getFinalTime() - pb.getInitialTime();
115          double scalAbsoluteTolerance = 1.0e-8;
116          double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
117    
118          FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
119                                                                     scalAbsoluteTolerance,
120                                                                     scalRelativeTolerance);
121          TestProblemHandler handler = new TestProblemHandler(pb, integ);
122          integ.addStepHandler(handler);
123          integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
124                          pb.getFinalTime(), new double[pb.getDimension()]);
125    
126          assertTrue(handler.getLastError() < 2.0e-7);
127          assertTrue(handler.getMaximalValueError() < 2.0e-7);
128          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
129          assertEquals("Dormand-Prince 5(4)", integ.getName());
130      }
131    
132      private static class DP54SmallLastHandler implements StepHandler {
133    
134        public DP54SmallLastHandler(double minStep) {
135          lastSeen = false;
136          this.minStep = minStep;
137        }
138    
139        public boolean requiresDenseOutput() {
140          return false;
141        }
142    
143        public void reset() {
144        }
145    
146        public void handleStep(StepInterpolator interpolator, boolean isLast) {
147          if (isLast) {
148            lastSeen = true;
149            double h = interpolator.getCurrentTime() - interpolator.getPreviousTime();
150            assertTrue(Math.abs(h) < minStep);
151          }
152        }
153    
154        public boolean wasLastSeen() {
155          return lastSeen;
156        }
157    
158        private boolean lastSeen;
159        private double  minStep;
160    
161      }
162    
163      public void testIncreasingTolerance()
164        throws DerivativeException, IntegratorException {
165    
166        int previousCalls = Integer.MAX_VALUE;
167        for (int i = -12; i < -2; ++i) {
168          TestProblem1 pb = new TestProblem1();
169          double minStep = 0;
170          double maxStep = pb.getFinalTime() - pb.getInitialTime();
171          double scalAbsoluteTolerance = Math.pow(10.0, i);
172          double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
173    
174          EmbeddedRungeKuttaIntegrator integ =
175              new DormandPrince54Integrator(minStep, maxStep,
176                                            scalAbsoluteTolerance, scalRelativeTolerance);
177          TestProblemHandler handler = new TestProblemHandler(pb, integ);
178          integ.setSafety(0.8);
179          integ.setMaxGrowth(5.0);
180          integ.setMinReduction(0.3);
181          integ.addStepHandler(handler);
182          integ.integrate(pb,
183                          pb.getInitialTime(), pb.getInitialState(),
184                          pb.getFinalTime(), new double[pb.getDimension()]);
185          assertEquals(0.8, integ.getSafety(), 1.0e-12);
186          assertEquals(5.0, integ.getMaxGrowth(), 1.0e-12);
187          assertEquals(0.3, integ.getMinReduction(), 1.0e-12);
188    
189          // the 0.7 factor is only valid for this test
190          // and has been obtained from trial and error
191          // there is no general relation between local and global errors
192          assertTrue(handler.getMaximalValueError() < (0.7 * scalAbsoluteTolerance));
193          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
194    
195          int calls = pb.getCalls();
196          assertEquals(integ.getEvaluations(), calls);
197          assertTrue(calls <= previousCalls);
198          previousCalls = calls;
199    
200        }
201    
202      }
203    
204      public void testEvents()
205        throws DerivativeException, IntegratorException {
206    
207        TestProblem4 pb = new TestProblem4();
208        double minStep = 0;
209        double maxStep = pb.getFinalTime() - pb.getInitialTime();
210        double scalAbsoluteTolerance = 1.0e-8;
211        double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
212    
213        FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
214                                                                   scalAbsoluteTolerance,
215                                                                   scalRelativeTolerance);
216        TestProblemHandler handler = new TestProblemHandler(pb, integ);
217        integ.addStepHandler(handler);
218        EventHandler[] functions = pb.getEventsHandlers();
219        for (int l = 0; l < functions.length; ++l) {
220          integ.addEventHandler(functions[l],
221                                     Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
222        }
223        assertEquals(functions.length, integ.getEventHandlers().size());
224        integ.integrate(pb,
225                        pb.getInitialTime(), pb.getInitialState(),
226                        pb.getFinalTime(), new double[pb.getDimension()]);
227    
228        assertTrue(handler.getMaximalValueError() < 5.0e-6);
229        assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
230        assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
231        integ.clearEventHandlers();
232        assertEquals(0, integ.getEventHandlers().size());
233    
234      }
235    
236      public void testKepler()
237        throws DerivativeException, IntegratorException {
238    
239        final TestProblem3 pb  = new TestProblem3(0.9);
240        double minStep = 0;
241        double maxStep = pb.getFinalTime() - pb.getInitialTime();
242        double scalAbsoluteTolerance = 1.0e-8;
243        double scalRelativeTolerance = scalAbsoluteTolerance;
244    
245        FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
246                                                                   scalAbsoluteTolerance,
247                                                                   scalRelativeTolerance);
248        integ.addStepHandler(new KeplerHandler(pb));
249        integ.integrate(pb,
250                        pb.getInitialTime(), pb.getInitialState(),
251                        pb.getFinalTime(), new double[pb.getDimension()]);
252    
253        assertEquals(integ.getEvaluations(), pb.getCalls());
254        assertTrue(pb.getCalls() < 2800);
255    
256      }
257    
258      public void testVariableSteps()
259        throws DerivativeException, IntegratorException {
260    
261        final TestProblem3 pb  = new TestProblem3(0.9);
262        double minStep = 0;
263        double maxStep = pb.getFinalTime() - pb.getInitialTime();
264        double scalAbsoluteTolerance = 1.0e-8;
265        double scalRelativeTolerance = scalAbsoluteTolerance;
266    
267        FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
268                                                                   scalAbsoluteTolerance,
269                                                                   scalRelativeTolerance);
270        integ.addStepHandler(new VariableHandler());
271        double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
272                                          pb.getFinalTime(), new double[pb.getDimension()]);
273        assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
274      }
275    
276      private static class KeplerHandler implements StepHandler {
277        public KeplerHandler(TestProblem3 pb) {
278          this.pb = pb;
279          reset();
280        }
281        public boolean requiresDenseOutput() {
282          return true;
283        }
284        public void reset() {
285          nbSteps = 0;
286          maxError = 0;
287        }
288        public void handleStep(StepInterpolator interpolator,
289                               boolean isLast)
290        throws DerivativeException {
291    
292          ++nbSteps;
293          for (int a = 1; a < 10; ++a) {
294    
295            double prev   = interpolator.getPreviousTime();
296            double curr   = interpolator.getCurrentTime();
297            double interp = ((10 - a) * prev + a * curr) / 10;
298            interpolator.setInterpolatedTime(interp);
299    
300            double[] interpolatedY = interpolator.getInterpolatedState ();
301            double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
302            double dx = interpolatedY[0] - theoreticalY[0];
303            double dy = interpolatedY[1] - theoreticalY[1];
304            double error = dx * dx + dy * dy;
305            if (error > maxError) {
306              maxError = error;
307            }
308          }
309          if (isLast) {
310            assertTrue(maxError < 7.0e-10);
311            assertTrue(nbSteps < 400);
312          }
313        }
314        private int nbSteps;
315        private double maxError;
316        private TestProblem3 pb;
317      }
318    
319      private static class VariableHandler implements StepHandler {
320        public VariableHandler() {
321          firstTime = true;
322          minStep = 0;
323          maxStep = 0;
324        }
325        public boolean requiresDenseOutput() {
326          return false;
327        }
328        public void reset() {
329          firstTime = true;
330          minStep = 0;
331          maxStep = 0;
332        }
333        public void handleStep(StepInterpolator interpolator,
334                               boolean isLast) {
335    
336          double step = Math.abs(interpolator.getCurrentTime()
337                                 - interpolator.getPreviousTime());
338          if (firstTime) {
339            minStep   = Math.abs(step);
340            maxStep   = minStep;
341            firstTime = false;
342          } else {
343            if (step < minStep) {
344              minStep = step;
345            }
346            if (step > maxStep) {
347              maxStep = step;
348            }
349          }
350    
351          if (isLast) {
352            assertTrue(minStep < (1.0 / 450.0));
353            assertTrue(maxStep > (1.0 / 4.2));
354          }
355        }  
356        private boolean firstTime;
357        private double  minStep;
358        private double  maxStep;
359      }
360    
361      public static Test suite() {
362        return new TestSuite(DormandPrince54IntegratorTest.class);
363      }
364    
365    }