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 org.apache.commons.math.ode.DerivativeException;
021 import org.apache.commons.math.ode.FirstOrderIntegrator;
022 import org.apache.commons.math.ode.IntegratorException;
023 import org.apache.commons.math.ode.TestProblem1;
024 import org.apache.commons.math.ode.TestProblem3;
025 import org.apache.commons.math.ode.TestProblem4;
026 import org.apache.commons.math.ode.TestProblem5;
027 import org.apache.commons.math.ode.TestProblemAbstract;
028 import org.apache.commons.math.ode.TestProblemHandler;
029 import org.apache.commons.math.ode.events.EventHandler;
030 import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
031 import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
032 import org.apache.commons.math.ode.sampling.StepHandler;
033 import org.apache.commons.math.ode.sampling.StepInterpolator;
034
035 import junit.framework.*;
036
037 public class DormandPrince54IntegratorTest
038 extends TestCase {
039
040 public DormandPrince54IntegratorTest(String name) {
041 super(name);
042 }
043
044 public void testDimensionCheck() {
045 try {
046 TestProblem1 pb = new TestProblem1();
047 DormandPrince54Integrator integrator = new DormandPrince54Integrator(0.0, 1.0,
048 1.0e-10, 1.0e-10);
049 integrator.integrate(pb,
050 0.0, new double[pb.getDimension()+10],
051 1.0, new double[pb.getDimension()+10]);
052 fail("an exception should have been thrown");
053 } catch(DerivativeException de) {
054 fail("wrong exception caught");
055 } catch(IntegratorException ie) {
056 }
057 }
058
059 public void testMinStep() {
060
061 try {
062 TestProblem1 pb = new TestProblem1();
063 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
064 double maxStep = pb.getFinalTime() - pb.getInitialTime();
065 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
066 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
067
068 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
069 vecAbsoluteTolerance,
070 vecRelativeTolerance);
071 TestProblemHandler handler = new TestProblemHandler(pb, integ);
072 integ.addStepHandler(handler);
073 integ.integrate(pb,
074 pb.getInitialTime(), pb.getInitialState(),
075 pb.getFinalTime(), new double[pb.getDimension()]);
076 fail("an exception should have been thrown");
077 } catch(DerivativeException de) {
078 fail("wrong exception caught");
079 } catch(IntegratorException ie) {
080 }
081
082 }
083
084 public void testSmallLastStep()
085 throws DerivativeException, IntegratorException {
086
087 TestProblemAbstract pb = new TestProblem5();
088 double minStep = 1.25;
089 double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime());
090 double scalAbsoluteTolerance = 6.0e-4;
091 double scalRelativeTolerance = 6.0e-4;
092
093 AdaptiveStepsizeIntegrator integ =
094 new DormandPrince54Integrator(minStep, maxStep,
095 scalAbsoluteTolerance,
096 scalRelativeTolerance);
097
098 DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
099 integ.addStepHandler(handler);
100 integ.setInitialStepSize(1.7);
101 integ.integrate(pb,
102 pb.getInitialTime(), pb.getInitialState(),
103 pb.getFinalTime(), new double[pb.getDimension()]);
104 assertTrue(handler.wasLastSeen());
105 assertEquals("Dormand-Prince 5(4)", integ.getName());
106
107 }
108
109 public void testBackward()
110 throws DerivativeException, IntegratorException {
111
112 TestProblem5 pb = new TestProblem5();
113 double minStep = 0;
114 double maxStep = pb.getFinalTime() - pb.getInitialTime();
115 double scalAbsoluteTolerance = 1.0e-8;
116 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
117
118 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
119 scalAbsoluteTolerance,
120 scalRelativeTolerance);
121 TestProblemHandler handler = new TestProblemHandler(pb, integ);
122 integ.addStepHandler(handler);
123 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
124 pb.getFinalTime(), new double[pb.getDimension()]);
125
126 assertTrue(handler.getLastError() < 2.0e-7);
127 assertTrue(handler.getMaximalValueError() < 2.0e-7);
128 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
129 assertEquals("Dormand-Prince 5(4)", integ.getName());
130 }
131
132 private static class DP54SmallLastHandler implements StepHandler {
133
134 public DP54SmallLastHandler(double minStep) {
135 lastSeen = false;
136 this.minStep = minStep;
137 }
138
139 public boolean requiresDenseOutput() {
140 return false;
141 }
142
143 public void reset() {
144 }
145
146 public void handleStep(StepInterpolator interpolator, boolean isLast) {
147 if (isLast) {
148 lastSeen = true;
149 double h = interpolator.getCurrentTime() - interpolator.getPreviousTime();
150 assertTrue(Math.abs(h) < minStep);
151 }
152 }
153
154 public boolean wasLastSeen() {
155 return lastSeen;
156 }
157
158 private boolean lastSeen;
159 private double minStep;
160
161 }
162
163 public void testIncreasingTolerance()
164 throws DerivativeException, IntegratorException {
165
166 int previousCalls = Integer.MAX_VALUE;
167 for (int i = -12; i < -2; ++i) {
168 TestProblem1 pb = new TestProblem1();
169 double minStep = 0;
170 double maxStep = pb.getFinalTime() - pb.getInitialTime();
171 double scalAbsoluteTolerance = Math.pow(10.0, i);
172 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
173
174 EmbeddedRungeKuttaIntegrator integ =
175 new DormandPrince54Integrator(minStep, maxStep,
176 scalAbsoluteTolerance, scalRelativeTolerance);
177 TestProblemHandler handler = new TestProblemHandler(pb, integ);
178 integ.setSafety(0.8);
179 integ.setMaxGrowth(5.0);
180 integ.setMinReduction(0.3);
181 integ.addStepHandler(handler);
182 integ.integrate(pb,
183 pb.getInitialTime(), pb.getInitialState(),
184 pb.getFinalTime(), new double[pb.getDimension()]);
185 assertEquals(0.8, integ.getSafety(), 1.0e-12);
186 assertEquals(5.0, integ.getMaxGrowth(), 1.0e-12);
187 assertEquals(0.3, integ.getMinReduction(), 1.0e-12);
188
189 // the 0.7 factor is only valid for this test
190 // and has been obtained from trial and error
191 // there is no general relation between local and global errors
192 assertTrue(handler.getMaximalValueError() < (0.7 * scalAbsoluteTolerance));
193 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
194
195 int calls = pb.getCalls();
196 assertEquals(integ.getEvaluations(), calls);
197 assertTrue(calls <= previousCalls);
198 previousCalls = calls;
199
200 }
201
202 }
203
204 public void testEvents()
205 throws DerivativeException, IntegratorException {
206
207 TestProblem4 pb = new TestProblem4();
208 double minStep = 0;
209 double maxStep = pb.getFinalTime() - pb.getInitialTime();
210 double scalAbsoluteTolerance = 1.0e-8;
211 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
212
213 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
214 scalAbsoluteTolerance,
215 scalRelativeTolerance);
216 TestProblemHandler handler = new TestProblemHandler(pb, integ);
217 integ.addStepHandler(handler);
218 EventHandler[] functions = pb.getEventsHandlers();
219 for (int l = 0; l < functions.length; ++l) {
220 integ.addEventHandler(functions[l],
221 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
222 }
223 assertEquals(functions.length, integ.getEventHandlers().size());
224 integ.integrate(pb,
225 pb.getInitialTime(), pb.getInitialState(),
226 pb.getFinalTime(), new double[pb.getDimension()]);
227
228 assertTrue(handler.getMaximalValueError() < 5.0e-6);
229 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
230 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
231 integ.clearEventHandlers();
232 assertEquals(0, integ.getEventHandlers().size());
233
234 }
235
236 public void testKepler()
237 throws DerivativeException, IntegratorException {
238
239 final TestProblem3 pb = new TestProblem3(0.9);
240 double minStep = 0;
241 double maxStep = pb.getFinalTime() - pb.getInitialTime();
242 double scalAbsoluteTolerance = 1.0e-8;
243 double scalRelativeTolerance = scalAbsoluteTolerance;
244
245 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
246 scalAbsoluteTolerance,
247 scalRelativeTolerance);
248 integ.addStepHandler(new KeplerHandler(pb));
249 integ.integrate(pb,
250 pb.getInitialTime(), pb.getInitialState(),
251 pb.getFinalTime(), new double[pb.getDimension()]);
252
253 assertEquals(integ.getEvaluations(), pb.getCalls());
254 assertTrue(pb.getCalls() < 2800);
255
256 }
257
258 public void testVariableSteps()
259 throws DerivativeException, IntegratorException {
260
261 final TestProblem3 pb = new TestProblem3(0.9);
262 double minStep = 0;
263 double maxStep = pb.getFinalTime() - pb.getInitialTime();
264 double scalAbsoluteTolerance = 1.0e-8;
265 double scalRelativeTolerance = scalAbsoluteTolerance;
266
267 FirstOrderIntegrator integ = new DormandPrince54Integrator(minStep, maxStep,
268 scalAbsoluteTolerance,
269 scalRelativeTolerance);
270 integ.addStepHandler(new VariableHandler());
271 double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
272 pb.getFinalTime(), new double[pb.getDimension()]);
273 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
274 }
275
276 private static class KeplerHandler implements StepHandler {
277 public KeplerHandler(TestProblem3 pb) {
278 this.pb = pb;
279 reset();
280 }
281 public boolean requiresDenseOutput() {
282 return true;
283 }
284 public void reset() {
285 nbSteps = 0;
286 maxError = 0;
287 }
288 public void handleStep(StepInterpolator interpolator,
289 boolean isLast)
290 throws DerivativeException {
291
292 ++nbSteps;
293 for (int a = 1; a < 10; ++a) {
294
295 double prev = interpolator.getPreviousTime();
296 double curr = interpolator.getCurrentTime();
297 double interp = ((10 - a) * prev + a * curr) / 10;
298 interpolator.setInterpolatedTime(interp);
299
300 double[] interpolatedY = interpolator.getInterpolatedState ();
301 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime());
302 double dx = interpolatedY[0] - theoreticalY[0];
303 double dy = interpolatedY[1] - theoreticalY[1];
304 double error = dx * dx + dy * dy;
305 if (error > maxError) {
306 maxError = error;
307 }
308 }
309 if (isLast) {
310 assertTrue(maxError < 7.0e-10);
311 assertTrue(nbSteps < 400);
312 }
313 }
314 private int nbSteps;
315 private double maxError;
316 private TestProblem3 pb;
317 }
318
319 private static class VariableHandler implements StepHandler {
320 public VariableHandler() {
321 firstTime = true;
322 minStep = 0;
323 maxStep = 0;
324 }
325 public boolean requiresDenseOutput() {
326 return false;
327 }
328 public void reset() {
329 firstTime = true;
330 minStep = 0;
331 maxStep = 0;
332 }
333 public void handleStep(StepInterpolator interpolator,
334 boolean isLast) {
335
336 double step = Math.abs(interpolator.getCurrentTime()
337 - interpolator.getPreviousTime());
338 if (firstTime) {
339 minStep = Math.abs(step);
340 maxStep = minStep;
341 firstTime = false;
342 } else {
343 if (step < minStep) {
344 minStep = step;
345 }
346 if (step > maxStep) {
347 maxStep = step;
348 }
349 }
350
351 if (isLast) {
352 assertTrue(minStep < (1.0 / 450.0));
353 assertTrue(maxStep > (1.0 / 4.2));
354 }
355 }
356 private boolean firstTime;
357 private double minStep;
358 private double maxStep;
359 }
360
361 public static Test suite() {
362 return new TestSuite(DormandPrince54IntegratorTest.class);
363 }
364
365 }