001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode.nonstiff;
019    
020    import junit.framework.Test;
021    import junit.framework.TestCase;
022    import junit.framework.TestSuite;
023    
024    import org.apache.commons.math.ConvergenceException;
025    import org.apache.commons.math.ode.DerivativeException;
026    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
027    import org.apache.commons.math.ode.FirstOrderIntegrator;
028    import org.apache.commons.math.ode.IntegratorException;
029    import org.apache.commons.math.ode.TestProblem1;
030    import org.apache.commons.math.ode.TestProblem3;
031    import org.apache.commons.math.ode.TestProblem4;
032    import org.apache.commons.math.ode.TestProblem5;
033    import org.apache.commons.math.ode.TestProblemHandler;
034    import org.apache.commons.math.ode.events.EventException;
035    import org.apache.commons.math.ode.events.EventHandler;
036    import org.apache.commons.math.ode.nonstiff.HighamHall54Integrator;
037    import org.apache.commons.math.ode.sampling.StepHandler;
038    import org.apache.commons.math.ode.sampling.StepInterpolator;
039    
040    public class HighamHall54IntegratorTest
041      extends TestCase {
042    
043      public HighamHall54IntegratorTest(String name) {
044        super(name);
045      }
046    
047      public void testWrongDerivative() {
048        try {
049          HighamHall54Integrator integrator =
050              new HighamHall54Integrator(0.0, 1.0, 1.0e-10, 1.0e-10);
051          FirstOrderDifferentialEquations equations =
052              new FirstOrderDifferentialEquations() {
053                private static final long serialVersionUID = -1157081786301178032L;
054                public void computeDerivatives(double t, double[] y, double[] dot)
055                throws DerivativeException {
056                if (t < -0.5) {
057                    throw new DerivativeException("{0}", "oops");
058                } else {
059                    throw new DerivativeException(new RuntimeException("oops"));
060               }
061              }
062              public int getDimension() {
063                  return 1;
064              }
065          };
066    
067          try  {
068            integrator.integrate(equations, -1.0, new double[1], 0.0, new double[1]);
069            fail("an exception should have been thrown");
070          } catch(DerivativeException de) {
071            // expected behavior
072          }
073    
074          try  {
075            integrator.integrate(equations, 0.0, new double[1], 1.0, new double[1]);
076            fail("an exception should have been thrown");
077          } catch(DerivativeException de) {
078            // expected behavior
079          }
080    
081        } catch (Exception e) {
082          fail("wrong exception caught: " + e.getMessage());        
083        }
084      }
085    
086      public void testMinStep() {
087    
088        try {
089          TestProblem1 pb = new TestProblem1();
090          double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
091          double maxStep = pb.getFinalTime() - pb.getInitialTime();
092          double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
093          double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
094    
095          FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
096                                                                  vecAbsoluteTolerance,
097                                                                  vecRelativeTolerance);
098          TestProblemHandler handler = new TestProblemHandler(pb, integ);
099          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    }