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 018 package org.apache.commons.math.ode.nonstiff; 019 020 import org.apache.commons.math.ode.DerivativeException; 021 import org.apache.commons.math.ode.FirstOrderDifferentialEquations; 022 import org.apache.commons.math.ode.IntegratorException; 023 import org.apache.commons.math.ode.events.CombinedEventsManager; 024 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 025 import org.apache.commons.math.ode.sampling.DummyStepInterpolator; 026 import org.apache.commons.math.ode.sampling.StepHandler; 027 028 /** 029 * This class implements the common part of all embedded Runge-Kutta 030 * integrators for Ordinary Differential Equations. 031 * 032 * <p>These methods are embedded explicit Runge-Kutta methods with two 033 * sets of coefficients allowing to estimate the error, their Butcher 034 * arrays are as follows : 035 * <pre> 036 * 0 | 037 * c2 | a21 038 * c3 | a31 a32 039 * ... | ... 040 * cs | as1 as2 ... ass-1 041 * |-------------------------- 042 * | b1 b2 ... bs-1 bs 043 * | b'1 b'2 ... b's-1 b's 044 * </pre> 045 * </p> 046 * 047 * <p>In fact, we rather use the array defined by ej = bj - b'j to 048 * compute directly the error rather than computing two estimates and 049 * then comparing them.</p> 050 * 051 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 052 * methods. This means the last evaluation of the derivatives in one 053 * step is the same as the first in the next step. Then, this 054 * evaluation can be reused from one step to the next one and the cost 055 * of such a method is really s-1 evaluations despite the method still 056 * has s stages. This behaviour is true only for successful steps, if 057 * the step is rejected after the error estimation phase, no 058 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 059 * asi = bi for all i.</p> 060 * 061 * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $ 062 * @since 1.2 063 */ 064 065 public abstract class EmbeddedRungeKuttaIntegrator 066 extends AdaptiveStepsizeIntegrator { 067 068 /** Indicator for <i>fsal</i> methods. */ 069 private final boolean fsal; 070 071 /** Time steps from Butcher array (without the first zero). */ 072 private final double[] c; 073 074 /** Internal weights from Butcher array (without the first empty row). */ 075 private final double[][] a; 076 077 /** External weights for the high order method from Butcher array. */ 078 private final double[] b; 079 080 /** Prototype of the step interpolator. */ 081 private final RungeKuttaStepInterpolator prototype; 082 083 /** Stepsize control exponent. */ 084 private final double exp; 085 086 /** Safety factor for stepsize control. */ 087 private double safety; 088 089 /** Minimal reduction factor for stepsize control. */ 090 private double minReduction; 091 092 /** Maximal growth factor for stepsize control. */ 093 private double maxGrowth; 094 095 /** Build a Runge-Kutta integrator with the given Butcher array. 096 * @param name name of the method 097 * @param fsal indicate that the method is an <i>fsal</i> 098 * @param c time steps from Butcher array (without the first zero) 099 * @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 210 // set up an interpolator sharing the integrator arrays 211 AbstractStepInterpolator interpolator; 212 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) { 213 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 214 rki.reinitialize(this, yTmp, yDotK, forward); 215 interpolator = rki; 216 } else { 217 interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward); 218 } 219 interpolator.storeTime(t0); 220 221 // set up integration control objects 222 stepStart = t0; 223 double hNew = 0; 224 boolean firstTime = true; 225 for (StepHandler handler : stepHandlers) { 226 handler.reset(); 227 } 228 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager); 229 boolean lastStep = false; 230 231 // main integration loop 232 while (!lastStep) { 233 234 interpolator.shift(); 235 236 double error = 0; 237 for (boolean loop = true; loop;) { 238 239 if (firstTime || !fsal) { 240 // first stage 241 computeDerivatives(stepStart, y, yDotK[0]); 242 } 243 244 if (firstTime) { 245 final double[] scale = new double[y0.length]; 246 if (vecAbsoluteTolerance == null) { 247 for (int i = 0; i < scale.length; ++i) { 248 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]); 249 } 250 } else { 251 for (int i = 0; i < scale.length; ++i) { 252 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * Math.abs(y[i]); 253 } 254 } 255 hNew = initializeStep(equations, forward, getOrder(), scale, 256 stepStart, y, yDotK[0], yTmp, yDotK[1]); 257 firstTime = false; 258 } 259 260 stepSize = hNew; 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 290 // discrete events handling 291 interpolator.storeTime(stepStart + stepSize); 292 if (manager.evaluateStep(interpolator)) { 293 final double dt = manager.getEventTime() - stepStart; 294 if (Math.abs(dt) <= Math.ulp(stepStart)) { 295 // we cannot simply truncate the step, reject the current computation 296 // and let the loop compute another state with the truncated step. 297 // it is so small (much probably exactly 0 due to limited accuracy) 298 // that the code above would fail handling it. 299 // So we set up an artificial 0 size step by copying states 300 interpolator.storeTime(stepStart); 301 System.arraycopy(y, 0, yTmp, 0, y0.length); 302 hNew = 0; 303 stepSize = 0; 304 loop = false; 305 } else { 306 // reject the step to match exactly the next switch time 307 hNew = dt; 308 } 309 } else { 310 // accept the step 311 loop = false; 312 } 313 314 } else { 315 // reject the step and attempt to reduce error by stepsize control 316 final double factor = 317 Math.min(maxGrowth, 318 Math.max(minReduction, safety * Math.pow(error, exp))); 319 hNew = filterStep(stepSize * factor, forward, false); 320 } 321 322 } 323 324 // the step has been accepted 325 final double nextStep = stepStart + stepSize; 326 System.arraycopy(yTmp, 0, y, 0, y0.length); 327 manager.stepAccepted(nextStep, y); 328 lastStep = manager.stop(); 329 330 // provide the step data to the step handler 331 interpolator.storeTime(nextStep); 332 for (StepHandler handler : stepHandlers) { 333 handler.handleStep(interpolator, lastStep); 334 } 335 stepStart = nextStep; 336 337 if (fsal) { 338 // save the last evaluation for the next step 339 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length); 340 } 341 342 if (manager.reset(stepStart, y) && ! lastStep) { 343 // some event handler has triggered changes that 344 // invalidate the derivatives, we need to recompute them 345 computeDerivatives(stepStart, y, yDotK[0]); 346 } 347 348 if (! lastStep) { 349 // in some rare cases we may get here with stepSize = 0, for example 350 // when an event occurs at integration start, reducing the first step 351 // to zero; we have to reset the step to some safe non zero value 352 stepSize = filterStep(stepSize, forward, true); 353 354 // stepsize control for next step 355 final double factor = Math.min(maxGrowth, 356 Math.max(minReduction, 357 safety * Math.pow(error, exp))); 358 final double scaledH = stepSize * factor; 359 final double nextT = stepStart + scaledH; 360 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 361 hNew = filterStep(scaledH, forward, nextIsLast); 362 } 363 364 } 365 366 final double stopTime = stepStart; 367 resetInternalState(); 368 return stopTime; 369 370 } 371 372 /** Get the minimal reduction factor for stepsize control. 373 * @return minimal reduction factor 374 */ 375 public double getMinReduction() { 376 return minReduction; 377 } 378 379 /** Set the minimal reduction factor for stepsize control. 380 * @param minReduction minimal reduction factor 381 */ 382 public void setMinReduction(final double minReduction) { 383 this.minReduction = minReduction; 384 } 385 386 /** Get the maximal growth factor for stepsize control. 387 * @return maximal growth factor 388 */ 389 public double getMaxGrowth() { 390 return maxGrowth; 391 } 392 393 /** Set the maximal growth factor for stepsize control. 394 * @param maxGrowth maximal growth factor 395 */ 396 public void setMaxGrowth(final double maxGrowth) { 397 this.maxGrowth = maxGrowth; 398 } 399 400 /** Compute the error ratio. 401 * @param yDotK derivatives computed during the first stages 402 * @param y0 estimate of the step at the start of the step 403 * @param y1 estimate of the step at the end of the step 404 * @param h current step 405 * @return error ratio, greater than 1 if step should be rejected 406 */ 407 protected abstract double estimateError(double[][] yDotK, 408 double[] y0, double[] y1, 409 double h); 410 411 }