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 018package org.apache.commons.math3.ode.nonstiff; 019 020import org.apache.commons.math3.exception.DimensionMismatchException; 021import org.apache.commons.math3.exception.MaxCountExceededException; 022import org.apache.commons.math3.exception.NoBracketingException; 023import org.apache.commons.math3.exception.NumberIsTooSmallException; 024import org.apache.commons.math3.ode.ExpandableStatefulODE; 025import org.apache.commons.math3.util.FastMath; 026 027/** 028 * This class implements the common part of all embedded Runge-Kutta 029 * integrators for Ordinary Differential Equations. 030 * 031 * <p>These methods are embedded explicit Runge-Kutta methods with two 032 * sets of coefficients allowing to estimate the error, their Butcher 033 * arrays are as follows : 034 * <pre> 035 * 0 | 036 * c2 | a21 037 * c3 | a31 a32 038 * ... | ... 039 * cs | as1 as2 ... ass-1 040 * |-------------------------- 041 * | b1 b2 ... bs-1 bs 042 * | b'1 b'2 ... b's-1 b's 043 * </pre> 044 * </p> 045 * 046 * <p>In fact, we rather use the array defined by ej = bj - b'j to 047 * compute directly the error rather than computing two estimates and 048 * then comparing them.</p> 049 * 050 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 051 * methods. This means the last evaluation of the derivatives in one 052 * step is the same as the first in the next step. Then, this 053 * evaluation can be reused from one step to the next one and the cost 054 * of such a method is really s-1 evaluations despite the method still 055 * has s stages. This behaviour is true only for successful steps, if 056 * the step is rejected after the error estimation phase, no 057 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 058 * asi = bi for all i.</p> 059 * 060 * @since 1.2 061 */ 062 063public abstract class EmbeddedRungeKuttaIntegrator 064 extends AdaptiveStepsizeIntegrator { 065 066 /** Indicator for <i>fsal</i> methods. */ 067 private final boolean fsal; 068 069 /** Time steps from Butcher array (without the first zero). */ 070 private final double[] c; 071 072 /** Internal weights from Butcher array (without the first empty row). */ 073 private final double[][] a; 074 075 /** External weights for the high order method from Butcher array. */ 076 private final double[] b; 077 078 /** Prototype of the step interpolator. */ 079 private final RungeKuttaStepInterpolator prototype; 080 081 /** Stepsize control exponent. */ 082 private final double exp; 083 084 /** Safety factor for stepsize control. */ 085 private double safety; 086 087 /** Minimal reduction factor for stepsize control. */ 088 private double minReduction; 089 090 /** Maximal growth factor for stepsize control. */ 091 private double maxGrowth; 092 093 /** Build a Runge-Kutta integrator with the given Butcher array. 094 * @param name name of the method 095 * @param fsal indicate that the method is an <i>fsal</i> 096 * @param c time steps from Butcher array (without the first zero) 097 * @param a internal weights from Butcher array (without the first empty row) 098 * @param b propagation weights for the high order method from Butcher array 099 * @param prototype prototype of the step interpolator to use 100 * @param minStep minimal step (sign is irrelevant, regardless of 101 * integration direction, forward or backward), the last step can 102 * be smaller than this 103 * @param maxStep maximal step (sign is irrelevant, regardless of 104 * integration direction, forward or backward), the last step can 105 * be smaller than this 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 void integrate(final ExpandableStatefulODE equations, final double t) 193 throws NumberIsTooSmallException, DimensionMismatchException, 194 MaxCountExceededException, NoBracketingException { 195 196 sanityChecks(equations, t); 197 setEquations(equations); 198 final boolean forward = t > equations.getTime(); 199 200 // create some internal working arrays 201 final double[] y0 = equations.getCompleteState(); 202 final double[] y = y0.clone(); 203 final int stages = c.length + 1; 204 final double[][] yDotK = new double[stages][y.length]; 205 final double[] yTmp = y0.clone(); 206 final double[] yDotTmp = new double[y.length]; 207 208 // set up an interpolator sharing the integrator arrays 209 final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy(); 210 interpolator.reinitialize(this, yTmp, yDotK, forward, 211 equations.getPrimaryMapper(), equations.getSecondaryMappers()); 212 interpolator.storeTime(equations.getTime()); 213 214 // set up integration control objects 215 stepStart = equations.getTime(); 216 double hNew = 0; 217 boolean firstTime = true; 218 initIntegration(equations.getTime(), y0, t); 219 220 // main integration loop 221 isLastStep = false; 222 do { 223 224 interpolator.shift(); 225 226 // iterate over step size, ensuring local normalized error is smaller than 1 227 double error = 10; 228 while (error >= 1.0) { 229 230 if (firstTime || !fsal) { 231 // first stage 232 computeDerivatives(stepStart, y, yDotK[0]); 233 } 234 235 if (firstTime) { 236 final double[] scale = new double[mainSetDimension]; 237 if (vecAbsoluteTolerance == null) { 238 for (int i = 0; i < scale.length; ++i) { 239 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]); 240 } 241 } else { 242 for (int i = 0; i < scale.length; ++i) { 243 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]); 244 } 245 } 246 hNew = initializeStep(forward, getOrder(), scale, 247 stepStart, y, yDotK[0], yTmp, yDotK[1]); 248 firstTime = false; 249 } 250 251 stepSize = hNew; 252 if (forward) { 253 if (stepStart + stepSize >= t) { 254 stepSize = t - stepStart; 255 } 256 } else { 257 if (stepStart + stepSize <= t) { 258 stepSize = t - stepStart; 259 } 260 } 261 262 // next stages 263 for (int k = 1; k < stages; ++k) { 264 265 for (int j = 0; j < y0.length; ++j) { 266 double sum = a[k-1][0] * yDotK[0][j]; 267 for (int l = 1; l < k; ++l) { 268 sum += a[k-1][l] * yDotK[l][j]; 269 } 270 yTmp[j] = y[j] + stepSize * sum; 271 } 272 273 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 274 275 } 276 277 // estimate the state at the end of the step 278 for (int j = 0; j < y0.length; ++j) { 279 double sum = b[0] * yDotK[0][j]; 280 for (int l = 1; l < stages; ++l) { 281 sum += b[l] * yDotK[l][j]; 282 } 283 yTmp[j] = y[j] + stepSize * sum; 284 } 285 286 // estimate the error at the end of the step 287 error = estimateError(yDotK, y, yTmp, stepSize); 288 if (error >= 1.0) { 289 // reject the step and attempt to reduce error by stepsize control 290 final double factor = 291 FastMath.min(maxGrowth, 292 FastMath.max(minReduction, safety * FastMath.pow(error, exp))); 293 hNew = filterStep(stepSize * factor, forward, false); 294 } 295 296 } 297 298 // local error is small enough: accept the step, trigger events and step handlers 299 interpolator.storeTime(stepStart + stepSize); 300 System.arraycopy(yTmp, 0, y, 0, y0.length); 301 System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); 302 stepStart = acceptStep(interpolator, y, yDotTmp, t); 303 System.arraycopy(y, 0, yTmp, 0, y.length); 304 305 if (!isLastStep) { 306 307 // prepare next step 308 interpolator.storeTime(stepStart); 309 310 if (fsal) { 311 // save the last evaluation for the next step 312 System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length); 313 } 314 315 // stepsize control for next step 316 final double factor = 317 FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp))); 318 final double scaledH = stepSize * factor; 319 final double nextT = stepStart + scaledH; 320 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 321 hNew = filterStep(scaledH, forward, nextIsLast); 322 323 final double filteredNextT = stepStart + hNew; 324 final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t); 325 if (filteredNextIsLast) { 326 hNew = t - stepStart; 327 } 328 329 } 330 331 } while (!isLastStep); 332 333 // dispatch results 334 equations.setTime(stepStart); 335 equations.setCompleteState(y); 336 337 resetInternalState(); 338 339 } 340 341 /** Get the minimal reduction factor for stepsize control. 342 * @return minimal reduction factor 343 */ 344 public double getMinReduction() { 345 return minReduction; 346 } 347 348 /** Set the minimal reduction factor for stepsize control. 349 * @param minReduction minimal reduction factor 350 */ 351 public void setMinReduction(final double minReduction) { 352 this.minReduction = minReduction; 353 } 354 355 /** Get the maximal growth factor for stepsize control. 356 * @return maximal growth factor 357 */ 358 public double getMaxGrowth() { 359 return maxGrowth; 360 } 361 362 /** Set the maximal growth factor for stepsize control. 363 * @param maxGrowth maximal growth factor 364 */ 365 public void setMaxGrowth(final double maxGrowth) { 366 this.maxGrowth = maxGrowth; 367 } 368 369 /** Compute the error ratio. 370 * @param yDotK derivatives computed during the first stages 371 * @param y0 estimate of the step at the start of the step 372 * @param y1 estimate of the step at the end of the step 373 * @param h current step 374 * @return error ratio, greater than 1 if step should be rejected 375 */ 376 protected abstract double estimateError(double[][] yDotK, 377 double[] y0, double[] y1, 378 double h); 379 380}