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.*;
21  
22  import org.apache.commons.math.ode.DerivativeException;
23  import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
24  import org.apache.commons.math.ode.FirstOrderIntegrator;
25  import org.apache.commons.math.ode.IntegratorException;
26  import org.apache.commons.math.ode.TestProblem1;
27  import org.apache.commons.math.ode.TestProblem3;
28  import org.apache.commons.math.ode.TestProblem5;
29  import org.apache.commons.math.ode.TestProblemAbstract;
30  import org.apache.commons.math.ode.TestProblemFactory;
31  import org.apache.commons.math.ode.TestProblemHandler;
32  import org.apache.commons.math.ode.events.EventHandler;
33  import org.apache.commons.math.ode.nonstiff.ThreeEighthesIntegrator;
34  import org.apache.commons.math.ode.sampling.StepHandler;
35  import org.apache.commons.math.ode.sampling.StepInterpolator;
36  
37  public class ThreeEighthesIntegratorTest
38    extends TestCase {
39  
40    public ThreeEighthesIntegratorTest(String name) {
41      super(name);
42    }
43  
44    public void testDimensionCheck() {
45      try  {
46        TestProblem1 pb = new TestProblem1();
47        new ThreeEighthesIntegrator(0.01).integrate(pb,
48                                                    0.0, new double[pb.getDimension()+10],
49                                                    1.0, new double[pb.getDimension()+10]);
50          fail("an exception should have been thrown");
51      } catch(DerivativeException de) {
52        fail("wrong exception caught");
53      } catch(IntegratorException ie) {
54      }
55    }
56    
57    public void testDecreasingSteps()
58      throws DerivativeException, IntegratorException  {
59        
60      TestProblemAbstract[] problems = TestProblemFactory.getProblems();
61      for (int k = 0; k < problems.length; ++k) {
62      
63        double previousError = Double.NaN;
64        for (int i = 4; i < 10; ++i) {
65  
66          TestProblemAbstract pb = problems[k].copy();
67          double step = (pb.getFinalTime() - pb.getInitialTime())
68            * Math.pow(2.0, -i);
69  
70          FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
71          TestProblemHandler handler = new TestProblemHandler(pb, integ);
72          integ.addStepHandler(handler);
73          EventHandler[] functions = pb.getEventsHandlers();
74          for (int l = 0; l < functions.length; ++l) {
75            integ.addEventHandler(functions[l],
76                                       Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
77          }
78          double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
79                                            pb.getFinalTime(), new double[pb.getDimension()]);
80          if (functions.length == 0) {
81              assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
82          }
83  
84          double error = handler.getMaximalValueError();
85          if (i > 4) {
86            assertTrue(error < Math.abs(previousError));
87          }
88          previousError = error;
89          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
90  
91        }
92  
93      }
94  
95    }
96  
97   public void testSmallStep()
98      throws DerivativeException, IntegratorException {
99  
100     TestProblem1 pb = new TestProblem1();
101     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
102 
103     FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
104     TestProblemHandler handler = new TestProblemHandler(pb, integ);
105     integ.addStepHandler(handler);
106     integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
107                     pb.getFinalTime(), new double[pb.getDimension()]);
108 
109     assertTrue(handler.getLastError() < 2.0e-13);
110     assertTrue(handler.getMaximalValueError() < 4.0e-12);
111     assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
112     assertEquals("3/8", integ.getName());
113 
114   }
115 
116   public void testBigStep()
117     throws DerivativeException, IntegratorException {
118 
119     TestProblem1 pb = new TestProblem1();
120     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
121 
122     FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
123     TestProblemHandler handler = new TestProblemHandler(pb, integ);
124     integ.addStepHandler(handler);
125     integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
126                     pb.getFinalTime(), new double[pb.getDimension()]);
127 
128     assertTrue(handler.getLastError() > 0.0004);
129     assertTrue(handler.getMaximalValueError() > 0.005);
130     assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
131 
132   }
133 
134   public void testBackward()
135       throws DerivativeException, IntegratorException {
136 
137       TestProblem5 pb = new TestProblem5();
138       double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
139 
140       FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
141       TestProblemHandler handler = new TestProblemHandler(pb, integ);
142       integ.addStepHandler(handler);
143       integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
144                       pb.getFinalTime(), new double[pb.getDimension()]);
145 
146       assertTrue(handler.getLastError() < 5.0e-10);
147       assertTrue(handler.getMaximalValueError() < 7.0e-10);
148       assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
149       assertEquals("3/8", integ.getName());
150   }
151 
152   public void testKepler()
153     throws DerivativeException, IntegratorException {
154 
155     final TestProblem3 pb  = new TestProblem3(0.9);
156     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
157 
158     FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
159     integ.addStepHandler(new KeplerHandler(pb));
160     integ.integrate(pb,
161                     pb.getInitialTime(), pb.getInitialState(),
162                     pb.getFinalTime(), new double[pb.getDimension()]);
163   }
164 
165   private static class KeplerHandler implements StepHandler {
166 
167     public KeplerHandler(TestProblem3 pb) {
168       this.pb = pb;
169       maxError = 0;
170     }
171 
172     public boolean requiresDenseOutput() {
173       return false;
174     }
175 
176     public void reset() {
177       maxError = 0;
178     }
179 
180     public void handleStep(StepInterpolator interpolator,
181                            boolean isLast) throws DerivativeException {
182 
183       double[] interpolatedY = interpolator.getInterpolatedState();
184       double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getCurrentTime());
185       double dx = interpolatedY[0] - theoreticalY[0];
186       double dy = interpolatedY[1] - theoreticalY[1];
187       double error = dx * dx + dy * dy;
188       if (error > maxError) {
189         maxError = error;
190       }
191       if (isLast) {
192         // even with more than 1000 evaluations per period,
193         // RK4 is not able to integrate such an eccentric
194         // orbit with a good accuracy
195         assertTrue(maxError > 0.005);
196       }
197     }
198 
199     private TestProblem3 pb;
200     private double maxError = 0;
201 
202   }
203 
204   public void testStepSize()
205     throws DerivativeException, IntegratorException {
206       final double step = 1.23456;
207       FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
208       integ.addStepHandler(new StepHandler() {
209           public void handleStep(StepInterpolator interpolator, boolean isLast) {
210               if (! isLast) {
211                   assertEquals(step,
212                                interpolator.getCurrentTime() - interpolator.getPreviousTime(),
213                                1.0e-12);
214               }
215           }
216           public boolean requiresDenseOutput() {
217               return false;
218           }
219           public void reset() {
220           }          
221       });
222       integ.integrate(new FirstOrderDifferentialEquations() {
223           public void computeDerivatives(double t, double[] y, double[] dot) {
224               dot[0] = 1.0;
225           }
226           public int getDimension() {
227               return 1;
228           }
229       }, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
230   }
231 
232   public static Test suite() {
233     return new TestSuite(ThreeEighthesIntegratorTest.class);
234   }
235 
236 }