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