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.*;
021    
022    import org.apache.commons.math.ode.DerivativeException;
023    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
024    import org.apache.commons.math.ode.FirstOrderIntegrator;
025    import org.apache.commons.math.ode.IntegratorException;
026    import org.apache.commons.math.ode.TestProblem1;
027    import org.apache.commons.math.ode.TestProblem3;
028    import org.apache.commons.math.ode.TestProblem5;
029    import org.apache.commons.math.ode.TestProblemAbstract;
030    import org.apache.commons.math.ode.TestProblemFactory;
031    import org.apache.commons.math.ode.TestProblemHandler;
032    import org.apache.commons.math.ode.events.EventHandler;
033    import org.apache.commons.math.ode.nonstiff.ClassicalRungeKuttaIntegrator;
034    import org.apache.commons.math.ode.sampling.StepHandler;
035    import org.apache.commons.math.ode.sampling.StepInterpolator;
036    
037    public class ClassicalRungeKuttaIntegratorTest
038      extends TestCase {
039    
040      public ClassicalRungeKuttaIntegratorTest(String name) {
041        super(name);
042      }
043    
044      public void testSanityChecks() {
045        try  {
046          TestProblem1 pb = new TestProblem1();
047          new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
048                                                            0.0, new double[pb.getDimension()+10],
049                                                            1.0, new double[pb.getDimension()]);
050            fail("an exception should have been thrown");
051        } catch(DerivativeException de) {
052          fail("wrong exception caught");
053        } catch(IntegratorException ie) {
054        }
055        try  {
056            TestProblem1 pb = new TestProblem1();
057            new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
058                                                              0.0, new double[pb.getDimension()],
059                                                              1.0, new double[pb.getDimension()+10]);
060              fail("an exception should have been thrown");
061          } catch(DerivativeException de) {
062            fail("wrong exception caught");
063          } catch(IntegratorException ie) {
064          }
065        try  {
066          TestProblem1 pb = new TestProblem1();
067          new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
068                                                            0.0, new double[pb.getDimension()],
069                                                            0.0, new double[pb.getDimension()]);
070            fail("an exception should have been thrown");
071        } catch(DerivativeException de) {
072          fail("wrong exception caught");
073        } catch(IntegratorException ie) {
074        }
075      }
076      
077      public void testDecreasingSteps()
078        throws DerivativeException, IntegratorException  {
079          
080        TestProblemAbstract[] problems = TestProblemFactory.getProblems();
081        for (int k = 0; k < problems.length; ++k) {
082    
083          double previousError = Double.NaN;
084          for (int i = 4; i < 10; ++i) {
085    
086            TestProblemAbstract pb = problems[k].copy();
087            double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i);
088    
089            FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
090            TestProblemHandler handler = new TestProblemHandler(pb, integ);
091            integ.addStepHandler(handler);
092            EventHandler[] functions = pb.getEventsHandlers();
093            for (int l = 0; l < functions.length; ++l) {
094              integ.addEventHandler(functions[l],
095                                         Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
096            }
097            assertEquals(functions.length, integ.getEventHandlers().size());
098            double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
099                                              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    }