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