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