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