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.GillIntegrator;
34  import org.apache.commons.math.ode.sampling.StepHandler;
35  import org.apache.commons.math.ode.sampling.StepInterpolator;
36  
37  public class GillIntegratorTest
38    extends TestCase {
39  
40    public GillIntegratorTest(String name) {
41      super(name);
42    }
43  
44    public void testDimensionCheck() {
45      try  {
46        TestProblem1 pb = new TestProblem1();
47        new GillIntegrator(0.01).integrate(pb,
48                                           0.0, new double[pb.getDimension()+10],
49                                           1.0, new double[pb.getDimension()+10]);
50          fail("an exception should have been thrown");
51      } catch(DerivativeException de) {
52        fail("wrong exception caught");
53      } catch(IntegratorException ie) {
54      }
55    }
56    
57    public void testDecreasingSteps()
58      throws DerivativeException, IntegratorException  {
59        
60      TestProblemAbstract[] problems = TestProblemFactory.getProblems();
61      for (int k = 0; k < problems.length; ++k) {
62  
63        double previousError = Double.NaN;
64        for (int i = 5; i < 10; ++i) {
65  
66          TestProblemAbstract pb = problems[k].copy();
67          double step = (pb.getFinalTime() - pb.getInitialTime())
68            * Math.pow(2.0, -i);
69  
70          FirstOrderIntegrator integ = new GillIntegrator(step);
71          TestProblemHandler handler = new TestProblemHandler(pb, integ);
72          integ.addStepHandler(handler);
73          EventHandler[] functions = pb.getEventsHandlers();
74          for (int l = 0; l < functions.length; ++l) {
75            integ.addEventHandler(functions[l],
76                                       Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
77          }
78          double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
79                                            pb.getFinalTime(), new double[pb.getDimension()]);
80          if (functions.length == 0) {
81              assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
82          }
83  
84          double error = handler.getMaximalValueError();
85          if (i > 5) {
86            assertTrue(error < Math.abs(previousError));
87          }
88          previousError = error;
89          assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
90  
91        }
92  
93      }
94  
95    }
96  
97    public void testSmallStep()
98      throws DerivativeException, IntegratorException {
99  
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 }