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 java.util.Arrays; 021 022 import org.apache.commons.math.ode.DerivativeException; 023 import org.apache.commons.math.ode.FirstOrderDifferentialEquations; 024 import org.apache.commons.math.ode.IntegratorException; 025 import org.apache.commons.math.ode.events.CombinedEventsManager; 026 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 027 import org.apache.commons.math.ode.sampling.DummyStepInterpolator; 028 import org.apache.commons.math.ode.sampling.StepHandler; 029 030 /** 031 * This class implements the common part of all embedded Runge-Kutta 032 * integrators for Ordinary Differential Equations. 033 * 034 * <p>These methods are embedded explicit Runge-Kutta methods with two 035 * sets of coefficients allowing to estimate the error, their Butcher 036 * arrays are as follows : 037 * <pre> 038 * 0 | 039 * c2 | a21 040 * c3 | a31 a32 041 * ... | ... 042 * cs | as1 as2 ... ass-1 043 * |-------------------------- 044 * | b1 b2 ... bs-1 bs 045 * | b'1 b'2 ... b's-1 b's 046 * </pre> 047 * </p> 048 * 049 * <p>In fact, we rather use the array defined by ej = bj - b'j to 050 * compute directly the error rather than computing two estimates and 051 * then comparing them.</p> 052 * 053 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 054 * methods. This means the last evaluation of the derivatives in one 055 * step is the same as the first in the next step. Then, this 056 * evaluation can be reused from one step to the next one and the cost 057 * of such a method is really s-1 evaluations despite the method still 058 * has s stages. This behaviour is true only for successful steps, if 059 * the step is rejected after the error estimation phase, no 060 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 061 * asi = bi for all i.</p> 062 * 063 * @version $Revision: 786874 $ $Date: 2009-06-20 14:09:16 -0400 (Sat, 20 Jun 2009) $ 064 * @since 1.2 065 */ 066 067 public abstract class EmbeddedRungeKuttaIntegrator 068 extends AdaptiveStepsizeIntegrator { 069 070 /** Build a Runge-Kutta integrator with the given Butcher array. 071 * @param name name of the method 072 * @param fsal indicate that the method is an <i>fsal</i> 073 * @param c time steps from Butcher array (without the first zero) 074 * @param a internal weights from Butcher array (without the first empty row) 075 * @param b propagation weights for the high order method from Butcher array 076 * @param prototype prototype of the step interpolator to use 077 * @param minStep minimal step (must be positive even for backward 078 * integration), the last step can be smaller than this 079 * @param maxStep maximal step (must be positive even for backward 080 * integration) 081 * @param scalAbsoluteTolerance allowed absolute error 082 * @param scalRelativeTolerance allowed relative error 083 */ 084 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 085 final double[] c, final double[][] a, final double[] b, 086 final RungeKuttaStepInterpolator prototype, 087 final double minStep, final double maxStep, 088 final double scalAbsoluteTolerance, 089 final double scalRelativeTolerance) { 090 091 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 092 093 this.fsal = fsal; 094 this.c = c; 095 this.a = a; 096 this.b = b; 097 this.prototype = prototype; 098 099 exp = -1.0 / getOrder(); 100 101 // set the default values of the algorithm control parameters 102 setSafety(0.9); 103 setMinReduction(0.2); 104 setMaxGrowth(10.0); 105 106 } 107 108 /** Build a Runge-Kutta integrator with the given Butcher array. 109 * @param name name of the method 110 * @param fsal indicate that the method is an <i>fsal</i> 111 * @param c time steps from Butcher array (without the first zero) 112 * @param a internal weights from Butcher array (without the first empty row) 113 * @param b propagation weights for the high order method from Butcher array 114 * @param prototype prototype of the step interpolator to use 115 * @param minStep minimal step (must be positive even for backward 116 * integration), the last step can be smaller than this 117 * @param maxStep maximal step (must be positive even for backward 118 * integration) 119 * @param vecAbsoluteTolerance allowed absolute error 120 * @param vecRelativeTolerance allowed relative error 121 */ 122 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal, 123 final double[] c, final double[][] a, final double[] b, 124 final RungeKuttaStepInterpolator prototype, 125 final double minStep, final double maxStep, 126 final double[] vecAbsoluteTolerance, 127 final double[] vecRelativeTolerance) { 128 129 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 130 131 this.fsal = fsal; 132 this.c = c; 133 this.a = a; 134 this.b = b; 135 this.prototype = prototype; 136 137 exp = -1.0 / getOrder(); 138 139 // set the default values of the algorithm control parameters 140 setSafety(0.9); 141 setMinReduction(0.2); 142 setMaxGrowth(10.0); 143 144 } 145 146 /** Get the order of the method. 147 * @return order of the method 148 */ 149 public abstract int getOrder(); 150 151 /** Get the safety factor for stepsize control. 152 * @return safety factor 153 */ 154 public double getSafety() { 155 return safety; 156 } 157 158 /** Set the safety factor for stepsize control. 159 * @param safety safety factor 160 */ 161 public void setSafety(final double safety) { 162 this.safety = safety; 163 } 164 165 /** {@inheritDoc} */ 166 @Override 167 public double integrate(final FirstOrderDifferentialEquations equations, 168 final double t0, final double[] y0, 169 final double t, final double[] y) 170 throws DerivativeException, IntegratorException { 171 172 sanityChecks(equations, t0, y0, t, y); 173 setEquations(equations); 174 resetEvaluations(); 175 final boolean forward = (t > t0); 176 177 // create some internal working arrays 178 final int stages = c.length + 1; 179 if (y != y0) { 180 System.arraycopy(y0, 0, y, 0, y0.length); 181 } 182 final double[][] yDotK = new double[stages][y0.length]; 183 final double[] yTmp = new double[y0.length]; 184 185 // set up an interpolator sharing the integrator arrays 186 AbstractStepInterpolator interpolator; 187 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) { 188 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy(); 189 rki.reinitialize(this, yTmp, yDotK, forward); 190 interpolator = rki; 191 } else { 192 interpolator = new DummyStepInterpolator(yTmp, forward); 193 } 194 interpolator.storeTime(t0); 195 196 // set up integration control objects 197 stepStart = t0; 198 double hNew = 0; 199 boolean firstTime = true; 200 for (StepHandler handler : stepHandlers) { 201 handler.reset(); 202 } 203 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager); 204 boolean lastStep = false; 205 206 // main integration loop 207 while (!lastStep) { 208 209 interpolator.shift(); 210 211 double error = 0; 212 for (boolean loop = true; loop;) { 213 214 if (firstTime || !fsal) { 215 // first stage 216 computeDerivatives(stepStart, y, yDotK[0]); 217 } 218 219 if (firstTime) { 220 final double[] scale; 221 if (vecAbsoluteTolerance != null) { 222 scale = vecAbsoluteTolerance; 223 } else { 224 scale = new double[y0.length]; 225 Arrays.fill(scale, scalAbsoluteTolerance); 226 } 227 hNew = initializeStep(equations, forward, getOrder(), scale, 228 stepStart, y, yDotK[0], yTmp, yDotK[1]); 229 firstTime = false; 230 } 231 232 stepSize = hNew; 233 234 // next stages 235 for (int k = 1; k < stages; ++k) { 236 237 for (int j = 0; j < y0.length; ++j) { 238 double sum = a[k-1][0] * yDotK[0][j]; 239 for (int l = 1; l < k; ++l) { 240 sum += a[k-1][l] * yDotK[l][j]; 241 } 242 yTmp[j] = y[j] + stepSize * sum; 243 } 244 245 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]); 246 247 } 248 249 // estimate the state at the end of the step 250 for (int j = 0; j < y0.length; ++j) { 251 double sum = b[0] * yDotK[0][j]; 252 for (int l = 1; l < stages; ++l) { 253 sum += b[l] * yDotK[l][j]; 254 } 255 yTmp[j] = y[j] + stepSize * sum; 256 } 257 258 // estimate the error at the end of the step 259 error = estimateError(yDotK, y, yTmp, stepSize); 260 if (error <= 1.0) { 261 262 // discrete events handling 263 interpolator.storeTime(stepStart + stepSize); 264 if (manager.evaluateStep(interpolator)) { 265 final double dt = manager.getEventTime() - stepStart; 266 if (Math.abs(dt) <= Math.ulp(stepStart)) { 267 // rejecting the step would lead to a too small next step, we accept it 268 loop = false; 269 } else { 270 // reject the step to match exactly the next switch time 271 hNew = dt; 272 } 273 } else { 274 // accept the step 275 loop = false; 276 } 277 278 } else { 279 // reject the step and attempt to reduce error by stepsize control 280 final double factor = 281 Math.min(maxGrowth, 282 Math.max(minReduction, safety * Math.pow(error, exp))); 283 hNew = filterStep(stepSize * factor, forward, false); 284 } 285 286 } 287 288 // the step has been accepted 289 final double nextStep = stepStart + stepSize; 290 System.arraycopy(yTmp, 0, y, 0, y0.length); 291 manager.stepAccepted(nextStep, y); 292 lastStep = manager.stop(); 293 294 // provide the step data to the step handler 295 interpolator.storeTime(nextStep); 296 for (StepHandler handler : stepHandlers) { 297 handler.handleStep(interpolator, lastStep); 298 } 299 stepStart = nextStep; 300 301 if (fsal) { 302 // save the last evaluation for the next step 303 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length); 304 } 305 306 if (manager.reset(stepStart, y) && ! lastStep) { 307 // some event handler has triggered changes that 308 // invalidate the derivatives, we need to recompute them 309 computeDerivatives(stepStart, y, yDotK[0]); 310 } 311 312 if (! lastStep) { 313 // in some rare cases we may get here with stepSize = 0, for example 314 // when an event occurs at integration start, reducing the first step 315 // to zero; we have to reset the step to some safe non zero value 316 stepSize = filterStep(stepSize, forward, true); 317 318 // stepsize control for next step 319 final double factor = Math.min(maxGrowth, 320 Math.max(minReduction, 321 safety * Math.pow(error, exp))); 322 final double scaledH = stepSize * factor; 323 final double nextT = stepStart + scaledH; 324 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); 325 hNew = filterStep(scaledH, forward, nextIsLast); 326 } 327 328 } 329 330 final double stopTime = stepStart; 331 resetInternalState(); 332 return stopTime; 333 334 } 335 336 /** Get the minimal reduction factor for stepsize control. 337 * @return minimal reduction factor 338 */ 339 public double getMinReduction() { 340 return minReduction; 341 } 342 343 /** Set the minimal reduction factor for stepsize control. 344 * @param minReduction minimal reduction factor 345 */ 346 public void setMinReduction(final double minReduction) { 347 this.minReduction = minReduction; 348 } 349 350 /** Get the maximal growth factor for stepsize control. 351 * @return maximal growth factor 352 */ 353 public double getMaxGrowth() { 354 return maxGrowth; 355 } 356 357 /** Set the maximal growth factor for stepsize control. 358 * @param maxGrowth maximal growth factor 359 */ 360 public void setMaxGrowth(final double maxGrowth) { 361 this.maxGrowth = maxGrowth; 362 } 363 364 /** Compute the error ratio. 365 * @param yDotK derivatives computed during the first stages 366 * @param y0 estimate of the step at the start of the step 367 * @param y1 estimate of the step at the end of the step 368 * @param h current step 369 * @return error ratio, greater than 1 if step should be rejected 370 */ 371 protected abstract double estimateError(double[][] yDotK, 372 double[] y0, double[] y1, 373 double h); 374 375 /** Indicator for <i>fsal</i> methods. */ 376 private boolean fsal; 377 378 /** Time steps from Butcher array (without the first zero). */ 379 private double[] c; 380 381 /** Internal weights from Butcher array (without the first empty row). */ 382 private double[][] a; 383 384 /** External weights for the high order method from Butcher array. */ 385 private double[] b; 386 387 /** Prototype of the step interpolator. */ 388 private RungeKuttaStepInterpolator prototype; 389 390 /** Stepsize control exponent. */ 391 private double exp; 392 393 /** Safety factor for stepsize control. */ 394 private double safety; 395 396 /** Minimal reduction factor for stepsize control. */ 397 private double minReduction; 398 399 /** Maximal growth factor for stepsize control. */ 400 private double maxGrowth; 401 402 }