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 }