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 }