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 }