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;
19  
20  import org.apache.commons.math.ode.AdaptiveStepsizeIntegrator;
21  import org.apache.commons.math.ode.DerivativeException;
22  import org.apache.commons.math.ode.DormandPrince54Integrator;
23  import org.apache.commons.math.ode.FirstOrderIntegrator;
24  import org.apache.commons.math.ode.IntegratorException;
25  import org.apache.commons.math.ode.StepHandler;
26  import org.apache.commons.math.ode.StepInterpolator;
27  import org.apache.commons.math.ode.SwitchingFunction;
28  
29  import junit.framework.*;
30  
31  public class DormandPrince54IntegratorTest
32    extends TestCase {
33  
34    public DormandPrince54IntegratorTest(String name) {
35      super(name);
36    }
37  
38    public void testDimensionCheck() {
39      try  {
40        TestProblem1 pb = new TestProblem1();
41        DormandPrince54Integrator integrator = new DormandPrince54Integrator(0.0, 1.0,
42                                                                             1.0e-10, 1.0e-10);
43        integrator.integrate(pb,
44                             0.0, new double[pb.getDimension()+10],
45                             1.0, new double[pb.getDimension()+10]);
46        fail("an exception should have been thrown");
47      } catch(DerivativeException de) {
48        fail("wrong exception caught");
49      } catch(IntegratorException ie) {
50      }
51    }
52  
53    public void testMinStep()
54      throws DerivativeException, IntegratorException {
55  
56      try {
57        TestProblem1 pb = new TestProblem1();
58        double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
59        double maxStep = pb.getFinalTime() - pb.getInitialTime();
60        double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
61        double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
62  
63        FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
64                                                                   vecAbsoluteTolerance,
65                                                                   vecRelativeTolerance);
66        TestProblemHandler handler = new TestProblemHandler(pb, integ);
67        integ.setStepHandler(handler);
68        integ.integrate(pb,
69                        pb.getInitialTime(), pb.getInitialState(),
70                        pb.getFinalTime(), new double[pb.getDimension()]);
71        fail("an exception should have been thrown");
72      } catch(DerivativeException de) {
73        fail("wrong exception caught");
74      } catch(IntegratorException ie) {
75      }
76  
77    }
78  
79    public void testSmallLastStep()
80      throws DerivativeException, IntegratorException {
81  
82      TestProblemAbstract pb = new TestProblem5();
83      double minStep = 1.25;
84      double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
85      double scalAbsoluteTolerance = 6.0e-4;
86      double scalRelativeTolerance = 6.0e-4;
87  
88      AdaptiveStepsizeIntegrator integ =
89        new DormandPrince54Integrator(minStep, maxStep,
90                                      scalAbsoluteTolerance,
91                                      scalRelativeTolerance);
92  
93      DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
94      integ.setStepHandler(handler);
95      integ.setInitialStepSize(1.7);
96      integ.integrate(pb,
97                      pb.getInitialTime(), pb.getInitialState(),
98                      pb.getFinalTime(), new double[pb.getDimension()]);
99      assertTrue(handler.wasLastSeen());
100     assertEquals("Dormand-Prince 5(4)", integ.getName());
101 
102   }
103 
104   private static class DP54SmallLastHandler implements StepHandler {
105 
106     public DP54SmallLastHandler(double minStep) {
107       lastSeen = false;
108       this.minStep = minStep;
109     }
110 
111     public boolean requiresDenseOutput() {
112       return false;
113     }
114 
115     public void reset() {
116     }
117 
118     public void handleStep(StepInterpolator interpolator, boolean isLast) {
119       if (isLast) {
120         lastSeen = true;
121         double h = interpolator.getCurrentTime() - interpolator.getPreviousTime();
122         assertTrue(Math.abs(h) < minStep);
123       }
124     }
125 
126     public boolean wasLastSeen() {
127       return lastSeen;
128     }
129 
130     private boolean lastSeen;
131     private double  minStep;
132 
133   }
134 
135   public void testIncreasingTolerance()
136     throws DerivativeException, IntegratorException {
137 
138     int previousCalls = Integer.MAX_VALUE;
139     for (int i = -12; i < -2; ++i) {
140       TestProblem1 pb = new TestProblem1();
141       double minStep = 0;
142       double maxStep = pb.getFinalTime() - pb.getInitialTime();
143       double scalAbsoluteTolerance = Math.pow(10.0, i);
144       double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
145 
146       EmbeddedRungeKuttaIntegrator integ =
147           new DormandPrince54Integrator(minStep, maxStep,
148                                         scalAbsoluteTolerance, scalRelativeTolerance);
149       TestProblemHandler handler = new TestProblemHandler(pb, integ);
150       integ.setSafety(0.8);
151       integ.setMaxGrowth(5.0);
152       integ.setMinReduction(0.3);
153       integ.setStepHandler(handler);
154       integ.integrate(pb,
155                       pb.getInitialTime(), pb.getInitialState(),
156                       pb.getFinalTime(), new double[pb.getDimension()]);
157       assertEquals(0.8, integ.getSafety(), 1.0e-12);
158       assertEquals(5.0, integ.getMaxGrowth(), 1.0e-12);
159       assertEquals(0.3, integ.getMinReduction(), 1.0e-12);
160 
161       // the 0.7 factor is only valid for this test
162       // and has been obtained from trial and error
163       // there is no general relation between local and global errors
164       assertTrue(handler.getMaximalValueError() < (0.7 * scalAbsoluteTolerance));
165       assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
166 
167       int calls = pb.getCalls();
168       assertTrue(calls <= previousCalls);
169       previousCalls = calls;
170 
171     }
172 
173   }
174 
175   public void testSwitchingFunctions()
176     throws DerivativeException, IntegratorException {
177 
178     TestProblem4 pb = new TestProblem4();
179     double minStep = 0;
180     double maxStep = pb.getFinalTime() - pb.getInitialTime();
181     double scalAbsoluteTolerance = 1.0e-8;
182     double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
183 
184     FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
185                                                                scalAbsoluteTolerance,
186                                                                scalRelativeTolerance);
187     TestProblemHandler handler = new TestProblemHandler(pb, integ);
188     integ.setStepHandler(handler);
189     SwitchingFunction[] functions = pb.getSwitchingFunctions();
190     for (int l = 0; l < functions.length; ++l) {
191       integ.addSwitchingFunction(functions[l],
192                                  Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
193     }
194     integ.integrate(pb,
195                     pb.getInitialTime(), pb.getInitialState(),
196                     pb.getFinalTime(), new double[pb.getDimension()]);
197 
198     assertTrue(handler.getMaximalValueError() < 5.0e-6);
199     assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
200     assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
201 
202   }
203 
204   public void testKepler()
205     throws DerivativeException, IntegratorException {
206 
207     final TestProblem3 pb  = new TestProblem3(0.9);
208     double minStep = 0;
209     double maxStep = pb.getFinalTime() - pb.getInitialTime();
210     double scalAbsoluteTolerance = 1.0e-8;
211     double scalRelativeTolerance = scalAbsoluteTolerance;
212 
213     FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
214                                                                scalAbsoluteTolerance,
215                                                                scalRelativeTolerance);
216     integ.setStepHandler(new KeplerHandler(pb));
217     integ.integrate(pb,
218                     pb.getInitialTime(), pb.getInitialState(),
219                     pb.getFinalTime(), new double[pb.getDimension()]);
220 
221     assertTrue(pb.getCalls() < 2800);
222 
223   }
224 
225   public void testVariableSteps()
226     throws DerivativeException, IntegratorException {
227 
228     final TestProblem3 pb  = new TestProblem3(0.9);
229     double minStep = 0;
230     double maxStep = pb.getFinalTime() - pb.getInitialTime();
231     double scalAbsoluteTolerance = 1.0e-8;
232     double scalRelativeTolerance = scalAbsoluteTolerance;
233 
234     FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
235                                                                scalAbsoluteTolerance,
236                                                                scalRelativeTolerance);
237     integ.setStepHandler(new VariableHandler());
238     integ.integrate(pb,
239                     pb.getInitialTime(), pb.getInitialState(),
240                     pb.getFinalTime(), new double[pb.getDimension()]);
241   }
242 
243   private static class KeplerHandler implements StepHandler {
244     public KeplerHandler(TestProblem3 pb) {
245       this.pb = pb;
246       reset();
247     }
248     public boolean requiresDenseOutput() {
249       return true;
250     }
251     public void reset() {
252       nbSteps = 0;
253       maxError = 0;
254     }
255     public void handleStep(StepInterpolator interpolator,
256                            boolean isLast)
257     throws DerivativeException {
258 
259       ++nbSteps;
260       for (int a = 1; a < 10; ++a) {
261 
262         double prev   = interpolator.getPreviousTime();
263         double curr   = interpolator.getCurrentTime();
264         double interp = ((10 - a) * prev + a * curr) / 10;
265         interpolator.setInterpolatedTime(interp);
266 
267         double[] interpolatedY = interpolator.getInterpolatedState ();
268         double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
269         double dx = interpolatedY[0] - theoreticalY[0];
270         double dy = interpolatedY[1] - theoreticalY[1];
271         double error = dx * dx + dy * dy;
272         if (error > maxError) {
273           maxError = error;
274         }
275       }
276       if (isLast) {
277         assertTrue(maxError < 7.0e-10);
278         assertTrue(nbSteps < 400);
279       }
280     }
281     private int nbSteps;
282     private double maxError;
283     private TestProblem3 pb;
284   }
285 
286   private static class VariableHandler implements StepHandler {
287     public VariableHandler() {
288       firstTime = true;
289       minStep = 0;
290       maxStep = 0;
291     }
292     public boolean requiresDenseOutput() {
293       return false;
294     }
295     public void reset() {
296       firstTime = true;
297       minStep = 0;
298       maxStep = 0;
299     }
300     public void handleStep(StepInterpolator interpolator,
301                            boolean isLast) {
302 
303       double step = Math.abs(interpolator.getCurrentTime()
304                              - interpolator.getPreviousTime());
305       if (firstTime) {
306         minStep   = Math.abs(step);
307         maxStep   = minStep;
308         firstTime = false;
309       } else {
310         if (step < minStep) {
311           minStep = step;
312         }
313         if (step > maxStep) {
314           maxStep = step;
315         }
316       }
317 
318       if (isLast) {
319         assertTrue(minStep < (1.0 / 450.0));
320         assertTrue(maxStep > (1.0 / 4.2));
321       }
322     }  
323     private boolean firstTime;
324     private double  minStep;
325     private double  maxStep;
326   }
327 
328   public static Test suite() {
329     return new TestSuite(DormandPrince54IntegratorTest.class);
330   }
331 
332 }