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.GillIntegrator;
034    import org.apache.commons.math.ode.sampling.StepHandler;
035    import org.apache.commons.math.ode.sampling.StepInterpolator;
036    
037    public class GillIntegratorTest
038      extends TestCase {
039    
040      public GillIntegratorTest(String name) {
041        super(name);
042      }
043    
044      public void testDimensionCheck() {
045        try  {
046          TestProblem1 pb = new TestProblem1();
047          new GillIntegrator(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 = 5; 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 GillIntegrator(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 > 5) {
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 GillIntegrator(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("Gill", 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 GillIntegrator(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 GillIntegrator(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("Gill", 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 GillIntegrator(step);
159        integ.addStepHandler(new KeplerStepHandler(pb));
160        integ.integrate(pb,
161                        pb.getInitialTime(), pb.getInitialState(),
162                        pb.getFinalTime(), new double[pb.getDimension()]);
163      }
164    
165      public void testUnstableDerivative()
166      throws DerivativeException, IntegratorException {
167        final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
168        FirstOrderIntegrator integ = new GillIntegrator(0.3);
169        integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
170        double[] y = { Double.NaN };
171        integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
172        assertEquals(8.0, y[0], 1.0e-12);
173      }
174    
175      private static class KeplerStepHandler implements StepHandler {
176        public KeplerStepHandler(TestProblem3 pb) {
177          this.pb = pb;
178          reset();
179        }
180        public boolean requiresDenseOutput() {
181          return false;
182        }
183        public void reset() {
184          maxError = 0;
185        }
186        public void handleStep(StepInterpolator interpolator,
187                               boolean isLast) throws DerivativeException {
188    
189          double[] interpolatedY = interpolator.getInterpolatedState();
190          double[] theoreticalY  = pb.computeTheoreticalState(interpolator.getCurrentTime());
191          double dx = interpolatedY[0] - theoreticalY[0];
192          double dy = interpolatedY[1] - theoreticalY[1];
193          double error = dx * dx + dy * dy;
194          if (error > maxError) {
195            maxError = error;
196          }
197          if (isLast) {
198            // even with more than 1000 evaluations per period,
199            // RK4 is not able to integrate such an eccentric
200            // orbit with a good accuracy
201            assertTrue(maxError > 0.001);
202          }
203        }
204        private double maxError;
205        private TestProblem3 pb;
206      }
207    
208      public void testStepSize()
209        throws DerivativeException, IntegratorException {
210          final double step = 1.23456;
211          FirstOrderIntegrator integ = new GillIntegrator(step);
212          integ.addStepHandler(new StepHandler() {
213              public void handleStep(StepInterpolator interpolator, boolean isLast) {
214                  if (! isLast) {
215                      assertEquals(step,
216                                   interpolator.getCurrentTime() - interpolator.getPreviousTime(),
217                                   1.0e-12);
218                  }
219              }
220              public boolean requiresDenseOutput() {
221                  return false;
222              }
223              public void reset() {
224              }          
225          });
226          integ.integrate(new FirstOrderDifferentialEquations() {
227              private static final long serialVersionUID = 0L;
228              public void computeDerivatives(double t, double[] y, double[] dot) {
229                  dot[0] = 1.0;
230              }
231              public int getDimension() {
232                  return 1;
233              }
234          }, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
235      }
236    
237      public static Test suite() {
238        return new TestSuite(GillIntegratorTest.class);
239      }
240    
241    }