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.DormandPrince54Integrator;
31  import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
32  import org.apache.commons.math.ode.sampling.StepHandler;
33  import org.apache.commons.math.ode.sampling.StepInterpolator;
34  
35  import junit.framework.*;
36  
37  public class DormandPrince54IntegratorTest
38    extends TestCase {
39  
40    public DormandPrince54IntegratorTest(String name) {
41      super(name);
42    }
43  
44    public void testDimensionCheck() {
45      try  {
46        TestProblem1 pb = new TestProblem1();
47        DormandPrince54Integrator integrator = new DormandPrince54Integrator(0.0, 1.0,
48                                                                             1.0e-10, 1.0e-10);
49        integrator.integrate(pb,
50                             0.0, new double[pb.getDimension()+10],
51                             1.0, new double[pb.getDimension()+10]);
52        fail("an exception should have been thrown");
53      } catch(DerivativeException de) {
54        fail("wrong exception caught");
55      } catch(IntegratorException ie) {
56      }
57    }
58  
59    public void testMinStep() {
60  
61      try {
62        TestProblem1 pb = new TestProblem1();
63        double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
64        double maxStep = pb.getFinalTime() - pb.getInitialTime();
65        double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
66        double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
67  
68        FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
69                                                                   vecAbsoluteTolerance,
70                                                                   vecRelativeTolerance);
71        TestProblemHandler handler = new TestProblemHandler(pb, integ);
72        integ.addStepHandler(handler);
73        integ.integrate(pb,
74                        pb.getInitialTime(), pb.getInitialState(),
75                        pb.getFinalTime(), new double[pb.getDimension()]);
76        fail("an exception should have been thrown");
77      } catch(DerivativeException de) {
78        fail("wrong exception caught");
79      } catch(IntegratorException ie) {
80      }
81  
82    }
83  
84    public void testSmallLastStep()
85      throws DerivativeException, IntegratorException {
86  
87      TestProblemAbstract pb = new TestProblem5();
88      double minStep = 1.25;
89      double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
90      double scalAbsoluteTolerance = 6.0e-4;
91      double scalRelativeTolerance = 6.0e-4;
92  
93      AdaptiveStepsizeIntegrator integ =
94        new DormandPrince54Integrator(minStep, maxStep,
95                                      scalAbsoluteTolerance,
96                                      scalRelativeTolerance);
97  
98      DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
99      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 }