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.math.ode; 19 20 /** 21 * This class implements the common part of all embedded Runge-Kutta 22 * integrators for Ordinary Differential Equations. 23 * 24 * <p>These methods are embedded explicit Runge-Kutta methods with two 25 * sets of coefficients allowing to estimate the error, their Butcher 26 * arrays are as follows : 27 * <pre> 28 * 0 | 29 * c2 | a21 30 * c3 | a31 a32 31 * ... | ... 32 * cs | as1 as2 ... ass-1 33 * |-------------------------- 34 * | b1 b2 ... bs-1 bs 35 * | b'1 b'2 ... b's-1 b's 36 * </pre> 37 * </p> 38 * 39 * <p>In fact, we rather use the array defined by ej = bj - b'j to 40 * compute directly the error rather than computing two estimates and 41 * then comparing them.</p> 42 * 43 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 44 * methods. This means the last evaluation of the derivatives in one 45 * step is the same as the first in the next step. Then, this 46 * evaluation can be reused from one step to the next one and the cost 47 * of such a method is really s-1 evaluations despite the method still 48 * has s stages. This behaviour is true only for successful steps, if 49 * the step is rejected after the error estimation phase, no 50 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 51 * asi = bi for all i.</p> 52 * 53 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 54 * @since 1.2 55 */ 56 57 public abstract class EmbeddedRungeKuttaIntegrator 58 extends AdaptiveStepsizeIntegrator { 59 60 /** Build a Runge-Kutta integrator with the given Butcher array. 61 * @param fsal indicate that the method is an <i>fsal</i> 62 * @param c time steps from Butcher array (without the first zero) 63 * @param a internal weights from Butcher array (without the first empty row) 64 * @param b propagation weights for the high order method from Butcher array 65 * @param prototype prototype of the step interpolator to use 66 * @param minStep minimal step (must be positive even for backward 67 * integration), the last step can be smaller than this 68 * @param maxStep maximal step (must be positive even for backward 69 * integration) 70 * @param scalAbsoluteTolerance allowed absolute error 71 * @param scalRelativeTolerance allowed relative error 72 */ 73 protected EmbeddedRungeKuttaIntegrator(boolean fsal, 74 double[] c, double[][] a, double[] b, 75 RungeKuttaStepInterpolator prototype, 76 double minStep, double maxStep, 77 double scalAbsoluteTolerance, 78 double scalRelativeTolerance) { 79 80 super(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 81 82 this.fsal = fsal; 83 this.c = c; 84 this.a = a; 85 this.b = b; 86 this.prototype = prototype; 87 88 exp = -1.0 / getOrder(); 89 90 // set the default values of the algorithm control parameters 91 setSafety(0.9); 92 setMinReduction(0.2); 93 setMaxGrowth(10.0); 94 95 } 96 97 /** Build a Runge-Kutta integrator with the given Butcher array. 98 * @param fsal indicate that the method is an <i>fsal</i> 99 * @param c time steps from Butcher array (without the first zero) 100 * @param a internal weights from Butcher array (without the first empty row) 101 * @param b propagation weights for the high order method from Butcher array 102 * @param prototype prototype of the step interpolator to use 103 * @param minStep minimal step (must be positive even for backward 104 * integration), the last step can be smaller than this 105 * @param maxStep maximal step (must be positive even for backward 106 * integration) 107 * @param vecAbsoluteTolerance allowed absolute error 108 * @param vecRelativeTolerance allowed relative error 109 */ 110 protected EmbeddedRungeKuttaIntegrator(boolean fsal, 111 double[] c, double[][] a, double[] b, 112 RungeKuttaStepInterpolator prototype, 113 double minStep, double maxStep, 114 double[] vecAbsoluteTolerance, 115 double[] vecRelativeTolerance) { 116 117 super(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 118 119 this.fsal = fsal; 120 this.c = c; 121 this.a = a; 122 this.b = b; 123 this.prototype = prototype; 124 125 exp = -1.0 / getOrder(); 126 127 // set the default values of the algorithm control parameters 128 setSafety(0.9); 129 setMinReduction(0.2); 130 setMaxGrowth(10.0); 131 132 } 133 134 /** Get the name of the method. 135 * @return name of the method 136 */ 137 public abstract String getName(); 138 139 /** Get the order of the method. 140 * @return order of the method 141 */ 142 public abstract int getOrder(); 143 144 /** Get the safety factor for stepsize control. 145 * @return safety factor 146 */ 147 public double getSafety() { 148 return safety; 149 } 150 151 /** Set the safety factor for stepsize control. 152 * @param safety safety factor 153 */ 154 public void setSafety(double safety) { 155 this.safety = safety; 156 } 157 158 /** Integrate the differential equations up to the given time. 159 * <p>This method solves an Initial Value Problem (IVP).</p> 160 * <p>Since this method stores some internal state variables made 161 * available in its public interface during integration ({@link 162 * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p> 163 * @param equations differential equations to integrate 164 * @param t0 initial time 165 * @param y0 initial value of the state vector at t0 166 * @param t target time for the integration 167 * (can be set to a value smaller than <code>t0</code> for backward integration) 168 * @param y placeholder where to put the state vector at each successful 169 * step (and hence at the end of integration), can be the same object as y0 170 * @throws IntegratorException if the integrator cannot perform integration 171 * @throws DerivativeException this exception is propagated to the caller if 172 * the underlying user function triggers one 173 */ 174 public void integrate(FirstOrderDifferentialEquations equations, 175 double t0, double[] y0, 176 double t, double[] y) 177 throws DerivativeException, IntegratorException { 178 179 sanityChecks(equations, t0, y0, t, y); 180 boolean forward = (t > t0); 181 182 // create some internal working arrays 183 int stages = c.length + 1; 184 if (y != y0) { 185 System.arraycopy(y0, 0, y, 0, y0.length); 186 } 187 double[][] yDotK = new double[stages][]; 188 for (int i = 0; i < stages; ++i) { 189 yDotK [i] = new double[y0.length]; 190 } 191 double[] yTmp = new double[y0.length]; 192 193 // set up an interpolator sharing the integrator arrays 194 AbstractStepInterpolator interpolator; 195 if (handler.requiresDenseOutput() || (! switchesHandler.isEmpty())) { 196 RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 197 rki.reinitialize(equations, yTmp, yDotK, forward); 198 interpolator = rki; 199 } else { 200 interpolator = new DummyStepInterpolator(yTmp, forward); 201 } 202 interpolator.storeTime(t0); 203 204 stepStart = t0; 205 double hNew = 0; 206 boolean firstTime = true; 207 boolean lastStep; 208 handler.reset(); 209 do { 210 211 interpolator.shift(); 212 213 double error = 0; 214 for (boolean loop = true; loop;) { 215 216 if (firstTime || !fsal) { 217 // first stage 218 equations.computeDerivatives(stepStart, y, yDotK[0]); 219 } 220 221 if (firstTime) { 222 double[] scale; 223 if (vecAbsoluteTolerance != null) { 224 scale = vecAbsoluteTolerance; 225 } else { 226 scale = new double[y0.length]; 227 for (int i = 0; i < scale.length; ++i) { 228 scale[i] = scalAbsoluteTolerance; 229 } 230 } 231 hNew = initializeStep(equations, forward, getOrder(), scale, 232 stepStart, y, yDotK[0], yTmp, yDotK[1]); 233 firstTime = false; 234 } 235 236 stepSize = hNew; 237 238 // step adjustment near bounds 239 if ((forward && (stepStart + stepSize > t)) || 240 ((! forward) && (stepStart + stepSize < t))) { 241 stepSize = t - stepStart; 242 } 243 244 // next stages 245 for (int k = 1; k < stages; ++k) { 246 247 for (int j = 0; j < y0.length; ++j) { 248 double sum = a[k-1][0] * yDotK[0][j]; 249 for (int l = 1; l < k; ++l) { 250 sum += a[k-1][l] * yDotK[l][j]; 251 } 252 yTmp[j] = y[j] + stepSize * sum; 253 } 254 255 equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 256 257 } 258 259 // estimate the state at the end of the step 260 for (int j = 0; j < y0.length; ++j) { 261 double sum = b[0] * yDotK[0][j]; 262 for (int l = 1; l < stages; ++l) { 263 sum += b[l] * yDotK[l][j]; 264 } 265 yTmp[j] = y[j] + stepSize * sum; 266 } 267 268 // estimate the error at the end of the step 269 error = estimateError(yDotK, y, yTmp, stepSize); 270 if (error <= 1.0) { 271 272 // Switching functions handling 273 interpolator.storeTime(stepStart + stepSize); 274 if (switchesHandler.evaluateStep(interpolator)) { 275 // reject the step to match exactly the next switch time 276 hNew = switchesHandler.getEventTime() - stepStart; 277 } else { 278 // accept the step 279 loop = false; 280 } 281 282 } else { 283 // reject the step and attempt to reduce error by stepsize control 284 double factor = Math.min(maxGrowth, 285 Math.max(minReduction, 286 safety * Math.pow(error, exp))); 287 hNew = filterStep(stepSize * factor, false); 288 } 289 290 } 291 292 // the step has been accepted 293 double nextStep = stepStart + stepSize; 294 System.arraycopy(yTmp, 0, y, 0, y0.length); 295 switchesHandler.stepAccepted(nextStep, y); 296 if (switchesHandler.stop()) { 297 lastStep = true; 298 } else { 299 lastStep = forward ? (nextStep >= t) : (nextStep <= t); 300 } 301 302 // provide the step data to the step handler 303 interpolator.storeTime(nextStep); 304 handler.handleStep(interpolator, lastStep); 305 stepStart = nextStep; 306 307 if (fsal) { 308 // save the last evaluation for the next step 309 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length); 310 } 311 312 if (switchesHandler.reset(stepStart, y) && ! lastStep) { 313 // some switching function has triggered changes that 314 // invalidate the derivatives, we need to recompute them 315 equations.computeDerivatives(stepStart, y, yDotK[0]); 316 } 317 318 if (! lastStep) { 319 // stepsize control for next step 320 double factor = Math.min(maxGrowth, 321 Math.max(minReduction, 322 safety * Math.pow(error, exp))); 323 double scaledH = stepSize * factor; 324 double nextT = stepStart + scaledH; 325 boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 326 hNew = filterStep(scaledH, nextIsLast); 327 } 328 329 } while (! lastStep); 330 331 resetInternalState(); 332 333 } 334 335 /** Get the minimal reduction factor for stepsize control. 336 * @return minimal reduction factor 337 */ 338 public double getMinReduction() { 339 return minReduction; 340 } 341 342 /** Set the minimal reduction factor for stepsize control. 343 * @param minReduction minimal reduction factor 344 */ 345 public void setMinReduction(double minReduction) { 346 this.minReduction = minReduction; 347 } 348 349 /** Get the maximal growth factor for stepsize control. 350 * @return maximal growth factor 351 */ 352 public double getMaxGrowth() { 353 return maxGrowth; 354 } 355 356 /** Set the maximal growth factor for stepsize control. 357 * @param maxGrowth maximal growth factor 358 */ 359 public void setMaxGrowth(double maxGrowth) { 360 this.maxGrowth = maxGrowth; 361 } 362 363 /** Compute the error ratio. 364 * @param yDotK derivatives computed during the first stages 365 * @param y0 estimate of the step at the start of the step 366 * @param y1 estimate of the step at the end of the step 367 * @param h current step 368 * @return error ratio, greater than 1 if step should be rejected 369 */ 370 protected abstract double estimateError(double[][] yDotK, 371 double[] y0, double[] y1, 372 double h); 373 374 /** Indicator for <i>fsal</i> methods. */ 375 private boolean fsal; 376 377 /** Time steps from Butcher array (without the first zero). */ 378 private double[] c; 379 380 /** Internal weights from Butcher array (without the first empty row). */ 381 private double[][] a; 382 383 /** External weights for the high order method from Butcher array. */ 384 private double[] b; 385 386 /** Prototype of the step interpolator. */ 387 private RungeKuttaStepInterpolator prototype; 388 389 /** Stepsize control exponent. */ 390 private double exp; 391 392 /** Safety factor for stepsize control. */ 393 private double safety; 394 395 /** Minimal reduction factor for stepsize control. */ 396 private double minReduction; 397 398 /** Maximal growth factor for stepsize control. */ 399 private double maxGrowth; 400 401 }