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 }