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.math4.legacy.ode.nonstiff; 019 020import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 021import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 022import org.apache.commons.math4.legacy.exception.NoBracketingException; 023import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 024import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE; 025import org.apache.commons.math4.core.jdkmath.JdkMath; 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 * 045 * <p>In fact, we rather use the array defined by ej = bj - b'j to 046 * compute directly the error rather than computing two estimates and 047 * then comparing them.</p> 048 * 049 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 050 * methods. This means the last evaluation of the derivatives in one 051 * step is the same as the first in the next step. Then, this 052 * evaluation can be reused from one step to the next one and the cost 053 * of such a method is really s-1 evaluations despite the method still 054 * has s stages. This behaviour is true only for successful steps, if 055 * the step is rejected after the error estimation phase, no 056 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 057 * asi = bi for all i.</p> 058 * 059 * @since 1.2 060 */ 061 062public abstract class EmbeddedRungeKuttaIntegrator 063 extends AdaptiveStepsizeIntegrator { 064 065 /** Indicator for <i>fsal</i> methods. */ 066 private final boolean fsal; 067 068 /** Time steps from Butcher array (without the first zero). */ 069 private final double[] c; 070 071 /** Internal weights from Butcher array (without the first empty row). */ 072 private final double[][] a; 073 074 /** External weights for the high order method from Butcher array. */ 075 private final double[] b; 076 077 /** Prototype of the step interpolator. */ 078 private final RungeKuttaStepInterpolator prototype; 079 080 /** Stepsize control exponent. */ 081 private final double exp; 082 083 /** Safety factor for stepsize control. */ 084 private double safety; 085 086 /** Minimal reduction factor for stepsize control. */ 087 private double minReduction; 088 089 /** Maximal growth factor for stepsize control. */ 090 private double maxGrowth; 091 092 /** Build a Runge-Kutta integrator with the given Butcher array. 093 * @param name name of the method 094 * @param fsal indicate that the method is an <i>fsal</i> 095 * @param c time steps from Butcher array (without the first zero) 096 * @param a internal weights from Butcher array (without the first empty row) 097 * @param b propagation weights for the high order method from Butcher array 098 * @param prototype prototype of the step interpolator to use 099 * @param minStep minimal step (sign is irrelevant, regardless of 100 * integration direction, forward or backward), the last step can 101 * be smaller than this 102 * @param maxStep maximal step (sign is irrelevant, regardless of 103 * integration direction, forward or backward), the last step can 104 * be smaller than this 105 * @param scalAbsoluteTolerance allowed absolute error 106 * @param scalRelativeTolerance allowed relative error 107 */ 108 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 109 final double[] c, final double[][] a, final double[] b, 110 final RungeKuttaStepInterpolator prototype, 111 final double minStep, final double maxStep, 112 final double scalAbsoluteTolerance, 113 final double scalRelativeTolerance) { 114 115 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 116 117 this.fsal = fsal; 118 this.c = c; 119 this.a = a; 120 this.b = b; 121 this.prototype = prototype; 122 123 exp = -1.0 / getOrder(); 124 125 // set the default values of the algorithm control parameters 126 setSafety(0.9); 127 setMinReduction(0.2); 128 setMaxGrowth(10.0); 129 } 130 131 /** Build a Runge-Kutta integrator with the given Butcher array. 132 * @param name name of the method 133 * @param fsal indicate that the method is an <i>fsal</i> 134 * @param c time steps from Butcher array (without the first zero) 135 * @param a internal weights from Butcher array (without the first empty row) 136 * @param b propagation weights for the high order method from Butcher array 137 * @param prototype prototype of the step interpolator to use 138 * @param minStep minimal step (must be positive even for backward 139 * integration), the last step can be smaller than this 140 * @param maxStep maximal step (must be positive even for backward 141 * integration) 142 * @param vecAbsoluteTolerance allowed absolute error 143 * @param vecRelativeTolerance allowed relative error 144 */ 145 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 146 final double[] c, final double[][] a, final double[] b, 147 final RungeKuttaStepInterpolator prototype, 148 final double minStep, final double maxStep, 149 final double[] vecAbsoluteTolerance, 150 final double[] vecRelativeTolerance) { 151 152 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 153 154 this.fsal = fsal; 155 this.c = c; 156 this.a = a; 157 this.b = b; 158 this.prototype = prototype; 159 160 exp = -1.0 / getOrder(); 161 162 // set the default values of the algorithm control parameters 163 setSafety(0.9); 164 setMinReduction(0.2); 165 setMaxGrowth(10.0); 166 } 167 168 /** Get the order of the method. 169 * @return order of the method 170 */ 171 public abstract int getOrder(); 172 173 /** Get the safety factor for stepsize control. 174 * @return safety factor 175 */ 176 public double getSafety() { 177 return safety; 178 } 179 180 /** Set the safety factor for stepsize control. 181 * @param safety safety factor 182 */ 183 public void setSafety(final double safety) { 184 this.safety = safety; 185 } 186 187 /** {@inheritDoc} */ 188 @Override 189 public void integrate(final ExpandableStatefulODE equations, final double t) 190 throws NumberIsTooSmallException, DimensionMismatchException, 191 MaxCountExceededException, NoBracketingException { 192 193 sanityChecks(equations, t); 194 setEquations(equations); 195 final boolean forward = t > equations.getTime(); 196 197 // create some internal working arrays 198 final double[] y0 = equations.getCompleteState(); 199 final double[] y = y0.clone(); 200 final int stages = c.length + 1; 201 final double[][] yDotK = new double[stages][y.length]; 202 final double[] yTmp = y0.clone(); 203 final double[] yDotTmp = new double[y.length]; 204 205 // set up an interpolator sharing the integrator arrays 206 final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy(); 207 interpolator.reinitialize(this, yTmp, yDotK, forward, 208 equations.getPrimaryMapper(), equations.getSecondaryMappers()); 209 interpolator.storeTime(equations.getTime()); 210 211 // set up integration control objects 212 stepStart = equations.getTime(); 213 double hNew = 0; 214 boolean firstTime = true; 215 initIntegration(equations.getTime(), y0, t); 216 217 // main integration loop 218 isLastStep = false; 219 do { 220 221 interpolator.shift(); 222 223 // iterate over step size, ensuring local normalized error is smaller than 1 224 double error = 10; 225 while (error >= 1.0) { 226 227 if (firstTime || !fsal) { 228 // first stage 229 computeDerivatives(stepStart, y, yDotK[0]); 230 } 231 232 if (firstTime) { 233 final double[] scale = new double[mainSetDimension]; 234 if (vecAbsoluteTolerance == null) { 235 for (int i = 0; i < scale.length; ++i) { 236 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * JdkMath.abs(y[i]); 237 } 238 } else { 239 for (int i = 0; i < scale.length; ++i) { 240 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * JdkMath.abs(y[i]); 241 } 242 } 243 hNew = initializeStep(forward, getOrder(), scale, 244 stepStart, y, yDotK[0], yTmp, yDotK[1]); 245 firstTime = false; 246 } 247 248 stepSize = hNew; 249 if (forward) { 250 if (stepStart + stepSize >= t) { 251 stepSize = t - stepStart; 252 } 253 } else { 254 if (stepStart + stepSize <= t) { 255 stepSize = t - stepStart; 256 } 257 } 258 259 // next stages 260 for (int k = 1; k < stages; ++k) { 261 262 for (int j = 0; j < y0.length; ++j) { 263 double sum = a[k-1][0] * yDotK[0][j]; 264 for (int l = 1; l < k; ++l) { 265 sum += a[k-1][l] * yDotK[l][j]; 266 } 267 yTmp[j] = y[j] + stepSize * sum; 268 } 269 270 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 271 } 272 273 // estimate the state at the end of the step 274 for (int j = 0; j < y0.length; ++j) { 275 double sum = b[0] * yDotK[0][j]; 276 for (int l = 1; l < stages; ++l) { 277 sum += b[l] * yDotK[l][j]; 278 } 279 yTmp[j] = y[j] + stepSize * sum; 280 } 281 282 // estimate the error at the end of the step 283 error = estimateError(yDotK, y, yTmp, stepSize); 284 if (error >= 1.0) { 285 // reject the step and attempt to reduce error by stepsize control 286 final double factor = 287 JdkMath.min(maxGrowth, 288 JdkMath.max(minReduction, safety * JdkMath.pow(error, exp))); 289 hNew = filterStep(stepSize * factor, forward, false); 290 } 291 } 292 293 // local error is small enough: accept the step, trigger events and step handlers 294 interpolator.storeTime(stepStart + stepSize); 295 System.arraycopy(yTmp, 0, y, 0, y0.length); 296 System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); 297 stepStart = acceptStep(interpolator, y, yDotTmp, t); 298 System.arraycopy(y, 0, yTmp, 0, y.length); 299 300 if (!isLastStep) { 301 302 // prepare next step 303 interpolator.storeTime(stepStart); 304 305 if (fsal) { 306 // save the last evaluation for the next step 307 System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length); 308 } 309 310 // stepsize control for next step 311 final double factor = 312 JdkMath.min(maxGrowth, JdkMath.max(minReduction, safety * JdkMath.pow(error, exp))); 313 final double scaledH = stepSize * factor; 314 final double nextT = stepStart + scaledH; 315 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 316 hNew = filterStep(scaledH, forward, nextIsLast); 317 318 final double filteredNextT = stepStart + hNew; 319 final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t); 320 if (filteredNextIsLast) { 321 hNew = t - stepStart; 322 } 323 } 324 } while (!isLastStep); 325 326 // dispatch results 327 equations.setTime(stepStart); 328 equations.setCompleteState(y); 329 330 resetInternalState(); 331 } 332 333 /** Get the minimal reduction factor for stepsize control. 334 * @return minimal reduction factor 335 */ 336 public double getMinReduction() { 337 return minReduction; 338 } 339 340 /** Set the minimal reduction factor for stepsize control. 341 * @param minReduction minimal reduction factor 342 */ 343 public void setMinReduction(final double minReduction) { 344 this.minReduction = minReduction; 345 } 346 347 /** Get the maximal growth factor for stepsize control. 348 * @return maximal growth factor 349 */ 350 public double getMaxGrowth() { 351 return maxGrowth; 352 } 353 354 /** Set the maximal growth factor for stepsize control. 355 * @param maxGrowth maximal growth factor 356 */ 357 public void setMaxGrowth(final double maxGrowth) { 358 this.maxGrowth = maxGrowth; 359 } 360 361 /** Compute the error ratio. 362 * @param yDotK derivatives computed during the first stages 363 * @param y0 estimate of the step at the start of the step 364 * @param y1 estimate of the step at the end of the step 365 * @param h current step 366 * @return error ratio, greater than 1 if step should be rejected 367 */ 368 protected abstract double estimateError(double[][] yDotK, 369 double[] y0, double[] y1, 370 double h); 371}