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.math3.ode.nonstiff;
19 
20 
21 import org.apache.commons.math3.exception.DimensionMismatchException;
22 import org.apache.commons.math3.exception.MaxCountExceededException;
23 import org.apache.commons.math3.exception.NoBracketingException;
24 import org.apache.commons.math3.exception.NumberIsTooSmallException;
25 import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
26 import org.apache.commons.math3.ode.FirstOrderIntegrator;
27 import org.apache.commons.math3.ode.TestProblem1;
28 import org.apache.commons.math3.ode.TestProblem2;
29 import org.apache.commons.math3.ode.TestProblem3;
30 import org.apache.commons.math3.ode.TestProblem4;
31 import org.apache.commons.math3.ode.TestProblem5;
32 import org.apache.commons.math3.ode.TestProblem6;
33 import org.apache.commons.math3.ode.TestProblemAbstract;
34 import org.apache.commons.math3.ode.TestProblemHandler;
35 import org.apache.commons.math3.ode.events.EventHandler;
36 import org.apache.commons.math3.ode.sampling.StepHandler;
37 import org.apache.commons.math3.ode.sampling.StepInterpolator;
38 import org.apache.commons.math3.util.FastMath;
39 import org.junit.Assert;
40 import org.junit.Test;
41 
42 public class ClassicalRungeKuttaIntegratorTest {
43 
44   @Test
testMissedEndEvent()45   public void testMissedEndEvent()
46       throws DimensionMismatchException, NumberIsTooSmallException,
47              MaxCountExceededException, NoBracketingException {
48       final double   t0     = 1878250320.0000029;
49       final double   tEvent = 1878250379.9999986;
50       final double[] k      = { 1.0e-4, 1.0e-5, 1.0e-6 };
51       FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() {
52 
53           public int getDimension() {
54               return k.length;
55           }
56 
57           public void computeDerivatives(double t, double[] y, double[] yDot) {
58               for (int i = 0; i < y.length; ++i) {
59                   yDot[i] = k[i] * y[i];
60               }
61           }
62       };
63 
64       ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(60.0);
65 
66       double[] y0   = new double[k.length];
67       for (int i = 0; i < y0.length; ++i) {
68           y0[i] = i + 1;
69       }
70       double[] y    = new double[k.length];
71 
72       double finalT = integrator.integrate(ode, t0, y0, tEvent, y);
73       Assert.assertEquals(tEvent, finalT, 5.0e-6);
74       for (int i = 0; i < y.length; ++i) {
75           Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9);
76       }
77 
78       integrator.addEventHandler(new EventHandler() {
79 
80           public void init(double t0, double[] y0, double t) {
81           }
82 
83           public void resetState(double t, double[] y) {
84           }
85 
86           public double g(double t, double[] y) {
87               return t - tEvent;
88           }
89 
90           public Action eventOccurred(double t, double[] y, boolean increasing) {
91               Assert.assertEquals(tEvent, t, 5.0e-6);
92               return Action.CONTINUE;
93           }
94       }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
95       finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y);
96       Assert.assertEquals(tEvent + 120, finalT, 5.0e-6);
97       for (int i = 0; i < y.length; ++i) {
98           Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9);
99       }
100 
101   }
102 
103   @Test
testSanityChecks()104   public void testSanityChecks()
105       throws DimensionMismatchException, NumberIsTooSmallException,
106              MaxCountExceededException, NoBracketingException {
107     try  {
108       TestProblem1 pb = new TestProblem1();
109       new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
110                                                         0.0, new double[pb.getDimension()+10],
111                                                         1.0, new double[pb.getDimension()]);
112         Assert.fail("an exception should have been thrown");
113     } catch(DimensionMismatchException ie) {
114     }
115     try  {
116         TestProblem1 pb = new TestProblem1();
117         new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
118                                                           0.0, new double[pb.getDimension()],
119                                                           1.0, new double[pb.getDimension()+10]);
120           Assert.fail("an exception should have been thrown");
121       } catch(DimensionMismatchException ie) {
122       }
123     try  {
124       TestProblem1 pb = new TestProblem1();
125       new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
126                                                         0.0, new double[pb.getDimension()],
127                                                         0.0, new double[pb.getDimension()]);
128         Assert.fail("an exception should have been thrown");
129     } catch(NumberIsTooSmallException ie) {
130     }
131   }
132 
133   @Test
testDecreasingSteps()134   public void testDecreasingSteps()
135       throws DimensionMismatchException, NumberIsTooSmallException,
136              MaxCountExceededException, NoBracketingException {
137 
138     for (TestProblemAbstract pb : new TestProblemAbstract[] {
139         new TestProblem1(), new TestProblem2(), new TestProblem3(),
140         new TestProblem4(), new TestProblem5(), new TestProblem6()
141     }) {
142 
143       double previousValueError = Double.NaN;
144       double previousTimeError = Double.NaN;
145       for (int i = 4; i < 10; ++i) {
146 
147         double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
148 
149         FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
150         TestProblemHandler handler = new TestProblemHandler(pb, integ);
151         integ.addStepHandler(handler);
152         EventHandler[] functions = pb.getEventsHandlers();
153         for (int l = 0; l < functions.length; ++l) {
154           integ.addEventHandler(functions[l],
155                                      Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
156         }
157         Assert.assertEquals(functions.length, integ.getEventHandlers().size());
158         double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
159                                           pb.getFinalTime(), new double[pb.getDimension()]);
160         if (functions.length == 0) {
161             Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
162         }
163 
164         double error = handler.getMaximalValueError();
165         if (i > 4) {
166           Assert.assertTrue(error < 1.01 * FastMath.abs(previousValueError));
167         }
168         previousValueError = error;
169 
170         double timeError = handler.getMaximalTimeError();
171         if (i > 4) {
172           Assert.assertTrue(timeError <= FastMath.abs(previousTimeError));
173         }
174         previousTimeError = timeError;
175 
176         integ.clearEventHandlers();
177         Assert.assertEquals(0, integ.getEventHandlers().size());
178       }
179 
180     }
181 
182   }
183 
184   @Test
testSmallStep()185   public void testSmallStep()
186       throws DimensionMismatchException, NumberIsTooSmallException,
187              MaxCountExceededException, NoBracketingException {
188 
189     TestProblem1 pb = new TestProblem1();
190     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
191 
192     FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
193     TestProblemHandler handler = new TestProblemHandler(pb, integ);
194     integ.addStepHandler(handler);
195     integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
196                     pb.getFinalTime(), new double[pb.getDimension()]);
197 
198     Assert.assertTrue(handler.getLastError() < 2.0e-13);
199     Assert.assertTrue(handler.getMaximalValueError() < 4.0e-12);
200     Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
201     Assert.assertEquals("classical Runge-Kutta", integ.getName());
202   }
203 
204   @Test
205   public void testBigStep()
206       throws DimensionMismatchException, NumberIsTooSmallException,
207              MaxCountExceededException, NoBracketingException {
208 
209     TestProblem1 pb = new TestProblem1();
210     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
211 
212     FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
213     TestProblemHandler handler = new TestProblemHandler(pb, integ);
214     integ.addStepHandler(handler);
215     integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
216                     pb.getFinalTime(), new double[pb.getDimension()]);
217 
218     Assert.assertTrue(handler.getLastError() > 0.0004);
219     Assert.assertTrue(handler.getMaximalValueError() > 0.005);
220     Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
221 
222   }
223 
224   @Test
testBackward()225   public void testBackward()
226       throws DimensionMismatchException, NumberIsTooSmallException,
227              MaxCountExceededException, NoBracketingException {
228 
229     TestProblem5 pb = new TestProblem5();
230     double step = FastMath.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
231 
232     FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
233     TestProblemHandler handler = new TestProblemHandler(pb, integ);
234     integ.addStepHandler(handler);
235     integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
236                     pb.getFinalTime(), new double[pb.getDimension()]);
237 
238     Assert.assertTrue(handler.getLastError() < 5.0e-10);
239     Assert.assertTrue(handler.getMaximalValueError() < 7.0e-10);
240     Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
241     Assert.assertEquals("classical Runge-Kutta", integ.getName());
242   }
243 
244   @Test
245   public void testKepler()
246       throws DimensionMismatchException, NumberIsTooSmallException,
247              MaxCountExceededException, NoBracketingException {
248 
249     final TestProblem3 pb  = new TestProblem3(0.9);
250     double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
251 
252     FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
253     integ.addStepHandler(new KeplerHandler(pb));
254     integ.integrate(pb,
255                     pb.getInitialTime(), pb.getInitialState(),
256                     pb.getFinalTime(), new double[pb.getDimension()]);
257   }
258 
259   private static class KeplerHandler implements StepHandler {
260     public KeplerHandler(TestProblem3 pb) {
261       this.pb = pb;
262       maxError = 0;
263     }
264     public void init(double t0, double[] y0, double t) {
265       maxError = 0;
266     }
267     public void handleStep(StepInterpolator interpolator, boolean isLast)
268         throws MaxCountExceededException {
269 
270       double[] interpolatedY = interpolator.getInterpolatedState ();
271       double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getCurrentTime());
272       double dx = interpolatedY[0] - theoreticalY[0];
273       double dy = interpolatedY[1] - theoreticalY[1];
274       double error = dx * dx + dy * dy;
275       if (error > maxError) {
276         maxError = error;
277       }
278       if (isLast) {
279         // even with more than 1000 evaluations per period,
280         // RK4 is not able to integrate such an eccentric
281         // orbit with a good accuracy
282         Assert.assertTrue(maxError > 0.005);
283       }
284     }
285     private double maxError = 0;
286     private TestProblem3 pb;
287   }
288 
289   @Test
testStepSize()290   public void testStepSize()
291       throws DimensionMismatchException, NumberIsTooSmallException,
292              MaxCountExceededException, NoBracketingException {
293       final double step = 1.23456;
294       FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
295       integ.addStepHandler(new StepHandler() {
296           public void handleStep(StepInterpolator interpolator, boolean isLast) {
297               if (! isLast) {
298                   Assert.assertEquals(step,
299                                interpolator.getCurrentTime() - interpolator.getPreviousTime(),
300                                1.0e-12);
301               }
302           }
303           public void init(double t0, double[] y0, double t) {
304           }
305       });
306       integ.integrate(new FirstOrderDifferentialEquations() {
307           public void computeDerivatives(double t, double[] y, double[] dot) {
308               dot[0] = 1.0;
309           }
310           public int getDimension() {
311               return 1;
312           }
313       }, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
314   }
315 
316   @Test
testTooLargeFirstStep()317   public void testTooLargeFirstStep() {
318 
319       RungeKuttaIntegrator integ = new ClassicalRungeKuttaIntegrator(0.5);
320       final double start = 0.0;
321       final double end   = 0.001;
322       FirstOrderDifferentialEquations equations = new FirstOrderDifferentialEquations() {
323 
324           public int getDimension() {
325               return 1;
326           }
327 
328           public void computeDerivatives(double t, double[] y, double[] yDot) {
329               Assert.assertTrue(t >= FastMath.nextAfter(start, Double.NEGATIVE_INFINITY));
330               Assert.assertTrue(t <= FastMath.nextAfter(end,   Double.POSITIVE_INFINITY));
331               yDot[0] = -100.0 * y[0];
332           }
333 
334       };
335 
336       integ.integrate(equations, start, new double[] { 1.0 }, end, new double[1]);
337 
338   }
339 
340 }
341