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