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