1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math3.ode.nonstiff; 19 20 21 import org.apache.commons.math3.exception.DimensionMismatchException; 22 import org.apache.commons.math3.exception.MaxCountExceededException; 23 import org.apache.commons.math3.exception.NoBracketingException; 24 import org.apache.commons.math3.exception.NumberIsTooSmallException; 25 import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; 26 import org.apache.commons.math3.ode.FirstOrderIntegrator; 27 import org.apache.commons.math3.ode.TestProblem1; 28 import org.apache.commons.math3.ode.TestProblem2; 29 import org.apache.commons.math3.ode.TestProblem3; 30 import org.apache.commons.math3.ode.TestProblem4; 31 import org.apache.commons.math3.ode.TestProblem5; 32 import org.apache.commons.math3.ode.TestProblem6; 33 import org.apache.commons.math3.ode.TestProblemAbstract; 34 import org.apache.commons.math3.ode.TestProblemHandler; 35 import org.apache.commons.math3.ode.events.EventHandler; 36 import org.apache.commons.math3.ode.sampling.StepHandler; 37 import org.apache.commons.math3.ode.sampling.StepInterpolator; 38 import org.apache.commons.math3.util.FastMath; 39 import org.junit.Assert; 40 import org.junit.Test; 41 42 public class ClassicalRungeKuttaIntegratorTest { 43 44 @Test testMissedEndEvent()45 public void testMissedEndEvent() 46 throws DimensionMismatchException, NumberIsTooSmallException, 47 MaxCountExceededException, NoBracketingException { 48 final double t0 = 1878250320.0000029; 49 final double tEvent = 1878250379.9999986; 50 final double[] k = { 1.0e-4, 1.0e-5, 1.0e-6 }; 51 FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { 52 53 public int getDimension() { 54 return k.length; 55 } 56 57 public void computeDerivatives(double t, double[] y, double[] yDot) { 58 for (int i = 0; i < y.length; ++i) { 59 yDot[i] = k[i] * y[i]; 60 } 61 } 62 }; 63 64 ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(60.0); 65 66 double[] y0 = new double[k.length]; 67 for (int i = 0; i < y0.length; ++i) { 68 y0[i] = i + 1; 69 } 70 double[] y = new double[k.length]; 71 72 double finalT = integrator.integrate(ode, t0, y0, tEvent, y); 73 Assert.assertEquals(tEvent, finalT, 5.0e-6); 74 for (int i = 0; i < y.length; ++i) { 75 Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9); 76 } 77 78 integrator.addEventHandler(new EventHandler() { 79 80 public void init(double t0, double[] y0, double t) { 81 } 82 83 public void resetState(double t, double[] y) { 84 } 85 86 public double g(double t, double[] y) { 87 return t - tEvent; 88 } 89 90 public Action eventOccurred(double t, double[] y, boolean increasing) { 91 Assert.assertEquals(tEvent, t, 5.0e-6); 92 return Action.CONTINUE; 93 } 94 }, Double.POSITIVE_INFINITY, 1.0e-20, 100); 95 finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y); 96 Assert.assertEquals(tEvent + 120, finalT, 5.0e-6); 97 for (int i = 0; i < y.length; ++i) { 98 Assert.assertEquals(y0[i] * FastMath.exp(k[i] * (finalT - t0)), y[i], 1.0e-9); 99 } 100 101 } 102 103 @Test testSanityChecks()104 public void testSanityChecks() 105 throws DimensionMismatchException, NumberIsTooSmallException, 106 MaxCountExceededException, NoBracketingException { 107 try { 108 TestProblem1 pb = new TestProblem1(); 109 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 110 0.0, new double[pb.getDimension()+10], 111 1.0, new double[pb.getDimension()]); 112 Assert.fail("an exception should have been thrown"); 113 } catch(DimensionMismatchException ie) { 114 } 115 try { 116 TestProblem1 pb = new TestProblem1(); 117 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 118 0.0, new double[pb.getDimension()], 119 1.0, new double[pb.getDimension()+10]); 120 Assert.fail("an exception should have been thrown"); 121 } catch(DimensionMismatchException ie) { 122 } 123 try { 124 TestProblem1 pb = new TestProblem1(); 125 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 126 0.0, new double[pb.getDimension()], 127 0.0, new double[pb.getDimension()]); 128 Assert.fail("an exception should have been thrown"); 129 } catch(NumberIsTooSmallException ie) { 130 } 131 } 132 133 @Test testDecreasingSteps()134 public void testDecreasingSteps() 135 throws DimensionMismatchException, NumberIsTooSmallException, 136 MaxCountExceededException, NoBracketingException { 137 138 for (TestProblemAbstract pb : new TestProblemAbstract[] { 139 new TestProblem1(), new TestProblem2(), new TestProblem3(), 140 new TestProblem4(), new TestProblem5(), new TestProblem6() 141 }) { 142 143 double previousValueError = Double.NaN; 144 double previousTimeError = Double.NaN; 145 for (int i = 4; i < 10; ++i) { 146 147 double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i); 148 149 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 150 TestProblemHandler handler = new TestProblemHandler(pb, integ); 151 integ.addStepHandler(handler); 152 EventHandler[] functions = pb.getEventsHandlers(); 153 for (int l = 0; l < functions.length; ++l) { 154 integ.addEventHandler(functions[l], 155 Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000); 156 } 157 Assert.assertEquals(functions.length, integ.getEventHandlers().size()); 158 double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 159 pb.getFinalTime(), new double[pb.getDimension()]); 160 if (functions.length == 0) { 161 Assert.assertEquals(pb.getFinalTime(), stopTime, 1.0e-10); 162 } 163 164 double error = handler.getMaximalValueError(); 165 if (i > 4) { 166 Assert.assertTrue(error < 1.01 * FastMath.abs(previousValueError)); 167 } 168 previousValueError = error; 169 170 double timeError = handler.getMaximalTimeError(); 171 if (i > 4) { 172 Assert.assertTrue(timeError <= FastMath.abs(previousTimeError)); 173 } 174 previousTimeError = timeError; 175 176 integ.clearEventHandlers(); 177 Assert.assertEquals(0, integ.getEventHandlers().size()); 178 } 179 180 } 181 182 } 183 184 @Test testSmallStep()185 public void testSmallStep() 186 throws DimensionMismatchException, NumberIsTooSmallException, 187 MaxCountExceededException, NoBracketingException { 188 189 TestProblem1 pb = new TestProblem1(); 190 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001; 191 192 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 193 TestProblemHandler handler = new TestProblemHandler(pb, integ); 194 integ.addStepHandler(handler); 195 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 196 pb.getFinalTime(), new double[pb.getDimension()]); 197 198 Assert.assertTrue(handler.getLastError() < 2.0e-13); 199 Assert.assertTrue(handler.getMaximalValueError() < 4.0e-12); 200 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 201 Assert.assertEquals("classical Runge-Kutta", integ.getName()); 202 } 203 204 @Test 205 public void testBigStep() 206 throws DimensionMismatchException, NumberIsTooSmallException, 207 MaxCountExceededException, NoBracketingException { 208 209 TestProblem1 pb = new TestProblem1(); 210 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2; 211 212 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 213 TestProblemHandler handler = new TestProblemHandler(pb, integ); 214 integ.addStepHandler(handler); 215 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 216 pb.getFinalTime(), new double[pb.getDimension()]); 217 218 Assert.assertTrue(handler.getLastError() > 0.0004); 219 Assert.assertTrue(handler.getMaximalValueError() > 0.005); 220 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 221 222 } 223 224 @Test testBackward()225 public void testBackward() 226 throws DimensionMismatchException, NumberIsTooSmallException, 227 MaxCountExceededException, NoBracketingException { 228 229 TestProblem5 pb = new TestProblem5(); 230 double step = FastMath.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001; 231 232 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 233 TestProblemHandler handler = new TestProblemHandler(pb, integ); 234 integ.addStepHandler(handler); 235 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 236 pb.getFinalTime(), new double[pb.getDimension()]); 237 238 Assert.assertTrue(handler.getLastError() < 5.0e-10); 239 Assert.assertTrue(handler.getMaximalValueError() < 7.0e-10); 240 Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 241 Assert.assertEquals("classical Runge-Kutta", integ.getName()); 242 } 243 244 @Test 245 public void testKepler() 246 throws DimensionMismatchException, NumberIsTooSmallException, 247 MaxCountExceededException, NoBracketingException { 248 249 final TestProblem3 pb = new TestProblem3(0.9); 250 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003; 251 252 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 253 integ.addStepHandler(new KeplerHandler(pb)); 254 integ.integrate(pb, 255 pb.getInitialTime(), pb.getInitialState(), 256 pb.getFinalTime(), new double[pb.getDimension()]); 257 } 258 259 private static class KeplerHandler implements StepHandler { 260 public KeplerHandler(TestProblem3 pb) { 261 this.pb = pb; 262 maxError = 0; 263 } 264 public void init(double t0, double[] y0, double t) { 265 maxError = 0; 266 } 267 public void handleStep(StepInterpolator interpolator, boolean isLast) 268 throws MaxCountExceededException { 269 270 double[] interpolatedY = interpolator.getInterpolatedState (); 271 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime()); 272 double dx = interpolatedY[0] - theoreticalY[0]; 273 double dy = interpolatedY[1] - theoreticalY[1]; 274 double error = dx * dx + dy * dy; 275 if (error > maxError) { 276 maxError = error; 277 } 278 if (isLast) { 279 // even with more than 1000 evaluations per period, 280 // RK4 is not able to integrate such an eccentric 281 // orbit with a good accuracy 282 Assert.assertTrue(maxError > 0.005); 283 } 284 } 285 private double maxError = 0; 286 private TestProblem3 pb; 287 } 288 289 @Test testStepSize()290 public void testStepSize() 291 throws DimensionMismatchException, NumberIsTooSmallException, 292 MaxCountExceededException, NoBracketingException { 293 final double step = 1.23456; 294 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 295 integ.addStepHandler(new StepHandler() { 296 public void handleStep(StepInterpolator interpolator, boolean isLast) { 297 if (! isLast) { 298 Assert.assertEquals(step, 299 interpolator.getCurrentTime() - interpolator.getPreviousTime(), 300 1.0e-12); 301 } 302 } 303 public void init(double t0, double[] y0, double t) { 304 } 305 }); 306 integ.integrate(new FirstOrderDifferentialEquations() { 307 public void computeDerivatives(double t, double[] y, double[] dot) { 308 dot[0] = 1.0; 309 } 310 public int getDimension() { 311 return 1; 312 } 313 }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); 314 } 315 316 @Test testTooLargeFirstStep()317 public void testTooLargeFirstStep() { 318 319 RungeKuttaIntegrator integ = new ClassicalRungeKuttaIntegrator(0.5); 320 final double start = 0.0; 321 final double end = 0.001; 322 FirstOrderDifferentialEquations equations = new FirstOrderDifferentialEquations() { 323 324 public int getDimension() { 325 return 1; 326 } 327 328 public void computeDerivatives(double t, double[] y, double[] yDot) { 329 Assert.assertTrue(t >= FastMath.nextAfter(start, Double.NEGATIVE_INFINITY)); 330 Assert.assertTrue(t <= FastMath.nextAfter(end, Double.POSITIVE_INFINITY)); 331 yDot[0] = -100.0 * y[0]; 332 } 333 334 }; 335 336 integ.integrate(equations, start, new double[] { 1.0 }, end, new double[1]); 337 338 } 339 340 } 341