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.ConvergenceException;
21  import org.apache.commons.math.FunctionEvaluationException;
22  import org.apache.commons.math.ode.DerivativeException;
23  import org.apache.commons.math.ode.FirstOrderIntegrator;
24  import org.apache.commons.math.ode.HighamHall54Integrator;
25  import org.apache.commons.math.ode.IntegratorException;
26  import org.apache.commons.math.ode.StepHandler;
27  import org.apache.commons.math.ode.StepInterpolator;
28  import org.apache.commons.math.ode.SwitchingFunction;
29  
30  import junit.framework.*;
31  
32  public class HighamHall54IntegratorTest
33    extends TestCase {
34  
35    public HighamHall54IntegratorTest(String name) {
36      super(name);
37    }
38  
39    public void testWrongDerivative() {
40      try {
41        HighamHall54Integrator integrator =
42            new HighamHall54Integrator(0.0, 1.0, 1.0e-10, 1.0e-10);
43        FirstOrderDifferentialEquations equations =
44            new FirstOrderDifferentialEquations() {
45            public void computeDerivatives(double t, double[] y, double[] dot)
46              throws DerivativeException {
47              if (t < -0.5) {
48                  throw new DerivativeException("{0}", new String[] { "oops" });
49              } else {
50                  throw new DerivativeException(new RuntimeException("oops"));
51             }
52            }
53            public int getDimension() {
54                return 1;
55            }
56        };
57  
58        try  {
59          integrator.integrate(equations, -1.0, new double[1], 0.0, new double[1]);
60          fail("an exception should have been thrown");
61        } catch(DerivativeException de) {
62          // expected behavior
63        }
64  
65        try  {
66          integrator.integrate(equations, 0.0, new double[1], 1.0, new double[1]);
67          fail("an exception should have been thrown");
68        } catch(DerivativeException de) {
69          // expected behavior
70        }
71  
72      } catch (Exception e) {
73        fail("wrong exception caught: " + e.getMessage());        
74      }
75    }
76  
77    public void testMinStep()
78      throws DerivativeException, IntegratorException {
79  
80      try {
81        TestProblem1 pb = new TestProblem1();
82        double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
83        double maxStep = pb.getFinalTime() - pb.getInitialTime();
84        double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
85        double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
86  
87        FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
88                                                                vecAbsoluteTolerance,
89                                                                vecRelativeTolerance);
90        TestProblemHandler handler = new TestProblemHandler(pb, integ);
91        integ.setStepHandler(handler);
92        integ.integrate(pb,
93                        pb.getInitialTime(), pb.getInitialState(),
94                        pb.getFinalTime(), new double[pb.getDimension()]);
95        fail("an exception should have been thrown");
96      } catch(DerivativeException de) {
97        fail("wrong exception caught");
98      } catch(IntegratorException ie) {
99      }
100 
101   }
102 
103   public void testIncreasingTolerance()
104     throws DerivativeException, IntegratorException {
105 
106     int previousCalls = Integer.MAX_VALUE;
107     for (int i = -12; i < -2; ++i) {
108       TestProblem1 pb = new TestProblem1();
109       double minStep = 0;
110       double maxStep = pb.getFinalTime() - pb.getInitialTime();
111       double scalAbsoluteTolerance = Math.pow(10.0, i);
112       double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
113 
114       FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
115                                                               scalAbsoluteTolerance,
116                                                               scalRelativeTolerance);
117       TestProblemHandler handler = new TestProblemHandler(pb, integ);
118       integ.setStepHandler(handler);
119       integ.integrate(pb,
120                       pb.getInitialTime(), pb.getInitialState(),
121                       pb.getFinalTime(), new double[pb.getDimension()]);
122 
123       // the 1.3 factor is only valid for this test
124       // and has been obtained from trial and error
125       // there is no general relation between local and global errors
126       assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
127       assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
128 
129       int calls = pb.getCalls();
130       assertTrue(calls <= previousCalls);
131       previousCalls = calls;
132 
133     }
134 
135   }
136 
137   public void testSwitchingFunctions()
138     throws DerivativeException, IntegratorException {
139 
140     TestProblem4 pb = new TestProblem4();
141     double minStep = 0;
142     double maxStep = pb.getFinalTime() - pb.getInitialTime();
143     double scalAbsoluteTolerance = 1.0e-8;
144     double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
145 
146     FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
147                                                             scalAbsoluteTolerance,
148                                                             scalRelativeTolerance);
149     TestProblemHandler handler = new TestProblemHandler(pb, integ);
150     integ.setStepHandler(handler);
151     SwitchingFunction[] functions = pb.getSwitchingFunctions();
152     for (int l = 0; l < functions.length; ++l) {
153       integ.addSwitchingFunction(functions[l],
154                                  Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
155     }
156     integ.integrate(pb,
157                     pb.getInitialTime(), pb.getInitialState(),
158                     pb.getFinalTime(), new double[pb.getDimension()]);
159 
160     assertTrue(handler.getMaximalValueError() < 1.0e-7);
161     assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
162     assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
163 
164   }
165 
166   public void testSwitchingFunctionsError()
167     throws DerivativeException, IntegratorException {
168 
169       final TestProblem1 pb = new TestProblem1();
170       double minStep = 0;
171       double maxStep = pb.getFinalTime() - pb.getInitialTime();
172       double scalAbsoluteTolerance = 1.0e-8;
173       double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
174 
175       FirstOrderIntegrator integ =
176           new HighamHall54Integrator(minStep, maxStep,
177                                      scalAbsoluteTolerance, scalRelativeTolerance);
178       TestProblemHandler handler = new TestProblemHandler(pb, integ);
179       integ.setStepHandler(handler);
180 
181       integ.addSwitchingFunction(new SwitchingFunction() {
182         public int eventOccurred(double t, double[] y) {
183           return SwitchingFunction.CONTINUE;
184         }
185         public double g(double t, double[] y) throws FunctionEvaluationException {
186           double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
187           double offset = t - middle;
188           if (offset > 0) {
189             throw new FunctionEvaluationException(t);
190           }
191           return offset;
192         }
193         public void resetState(double t, double[] y) {
194         }
195         private static final long serialVersionUID = 935652725339916361L;
196       }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
197 
198       try {
199         integ.integrate(pb,
200                         pb.getInitialTime(), pb.getInitialState(),
201                         pb.getFinalTime(), new double[pb.getDimension()]);
202         fail("an exception should have been thrown");
203       } catch (IntegratorException ie) {
204         // expected behavior
205       } catch (Exception e) {
206         fail("wrong exception type caught");
207       }
208 
209   }
210 
211   public void testSwitchingFunctionsNoConvergence()
212   throws DerivativeException, IntegratorException {
213 
214     final TestProblem1 pb = new TestProblem1();
215     double minStep = 0;
216     double maxStep = pb.getFinalTime() - pb.getInitialTime();
217     double scalAbsoluteTolerance = 1.0e-8;
218     double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
219 
220     FirstOrderIntegrator integ =
221         new HighamHall54Integrator(minStep, maxStep,
222                                    scalAbsoluteTolerance, scalRelativeTolerance);
223     TestProblemHandler handler = new TestProblemHandler(pb, integ);
224     integ.setStepHandler(handler);
225 
226     integ.addSwitchingFunction(new SwitchingFunction() {
227       public int eventOccurred(double t, double[] y) {
228         return SwitchingFunction.CONTINUE;
229       }
230       public double g(double t, double[] y) {
231         double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
232         double offset = t - middle;
233         return (offset > 0) ? (offset + 0.5) : (offset - 0.5);
234       }
235       public void resetState(double t, double[] y) {
236       }
237       private static final long serialVersionUID = 935652725339916361L;
238     }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
239 
240     try {
241       integ.integrate(pb,
242                       pb.getInitialTime(), pb.getInitialState(),
243                       pb.getFinalTime(), new double[pb.getDimension()]);
244       fail("an exception should have been thrown");
245     } catch (IntegratorException ie) {
246        assertTrue(ie.getCause() != null);
247        assertTrue(ie.getCause() instanceof ConvergenceException);
248     } catch (Exception e) {
249       fail("wrong exception type caught");
250     }
251 
252 }
253 
254   public void testSanityChecks() {
255     try {
256       final TestProblem3 pb  = new TestProblem3(0.9);
257       double minStep = 0;
258       double maxStep = pb.getFinalTime() - pb.getInitialTime();
259 
260       try {
261         FirstOrderIntegrator integ =
262             new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
263         integ.integrate(pb, pb.getInitialTime(), new double[6],
264                         pb.getFinalTime(), new double[pb.getDimension()]);
265         fail("an exception should have been thrown");
266       } catch (IntegratorException ie) {
267         // expected behavior
268       }
269 
270       try {
271         FirstOrderIntegrator integ =
272             new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
273         integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
274                         pb.getFinalTime(), new double[6]);
275         fail("an exception should have been thrown");
276       } catch (IntegratorException ie) {
277         // expected behavior
278       }
279 
280       try {
281         FirstOrderIntegrator integ =
282             new HighamHall54Integrator(minStep, maxStep, new double[2], new double[4]);
283         integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
284                         pb.getFinalTime(), new double[pb.getDimension()]);
285         fail("an exception should have been thrown");
286       } catch (IntegratorException ie) {
287         // expected behavior
288       }
289 
290       try {
291         FirstOrderIntegrator integ =
292             new HighamHall54Integrator(minStep, maxStep, new double[4], new double[2]);
293         integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
294                         pb.getFinalTime(), new double[pb.getDimension()]);
295         fail("an exception should have been thrown");
296       } catch (IntegratorException ie) {
297         // expected behavior
298       }
299 
300       try {
301         FirstOrderIntegrator integ =
302             new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
303         integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
304                         pb.getInitialTime(), new double[pb.getDimension()]);
305         fail("an exception should have been thrown");
306       } catch (IntegratorException ie) {
307         // expected behavior
308       }
309 
310     } catch (Exception e) {
311       fail("wrong exception caught: " + e.getMessage());
312     }
313   }
314 
315   public void testKepler()
316     throws DerivativeException, IntegratorException {
317 
318     final TestProblem3 pb  = new TestProblem3(0.9);
319     double minStep = 0;
320     double maxStep = pb.getFinalTime() - pb.getInitialTime();
321     double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
322     double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
323 
324     FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
325                                                             vecAbsoluteTolerance,
326                                                             vecRelativeTolerance);
327     integ.setStepHandler(new KeplerHandler(pb));
328     integ.integrate(pb,
329                     pb.getInitialTime(), pb.getInitialState(),
330                     pb.getFinalTime(), new double[pb.getDimension()]);
331     assertEquals("Higham-Hall 5(4)", integ.getName());
332   }
333 
334   private static class KeplerHandler implements StepHandler {
335     public KeplerHandler(TestProblem3 pb) {
336       this.pb = pb;
337       nbSteps = 0;
338       maxError = 0;
339     }
340     public boolean requiresDenseOutput() {
341       return false;
342     }
343     public void reset() {
344       nbSteps = 0;
345       maxError = 0;
346     }
347     public void handleStep(StepInterpolator interpolator,
348                            boolean isLast) {
349 
350       ++nbSteps;
351       double[] interpolatedY = interpolator.getInterpolatedState ();
352       double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getCurrentTime());
353       double dx = interpolatedY[0] - theoreticalY[0];
354       double dy = interpolatedY[1] - theoreticalY[1];
355       double error = dx * dx + dy * dy;
356       if (error > maxError) {
357         maxError = error;
358       }
359       if (isLast) {
360         assertTrue(maxError < 4e-11);
361         assertTrue(nbSteps < 670);
362       }
363     }
364     private TestProblem3 pb;
365     private int nbSteps;
366     private double maxError;
367   }
368 
369   public static Test suite() {
370     return new TestSuite(HighamHall54IntegratorTest.class);
371   }
372 
373 }