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