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