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