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.math3.ode.nonstiff;
019
020import org.apache.commons.math3.exception.DimensionMismatchException;
021import org.apache.commons.math3.exception.MaxCountExceededException;
022import org.apache.commons.math3.exception.NoBracketingException;
023import org.apache.commons.math3.exception.NumberIsTooSmallException;
024import org.apache.commons.math3.ode.ExpandableStatefulODE;
025import 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 * @since 1.2
061 */
062
063public abstract class EmbeddedRungeKuttaIntegrator
064  extends AdaptiveStepsizeIntegrator {
065
066    /** Indicator for <i>fsal</i> methods. */
067    private final boolean fsal;
068
069    /** Time steps from Butcher array (without the first zero). */
070    private final double[] c;
071
072    /** Internal weights from Butcher array (without the first empty row). */
073    private final double[][] a;
074
075    /** External weights for the high order method from Butcher array. */
076    private final double[] b;
077
078    /** Prototype of the step interpolator. */
079    private final RungeKuttaStepInterpolator prototype;
080
081    /** Stepsize control exponent. */
082    private final double exp;
083
084    /** Safety factor for stepsize control. */
085    private double safety;
086
087    /** Minimal reduction factor for stepsize control. */
088    private double minReduction;
089
090    /** Maximal growth factor for stepsize control. */
091    private double maxGrowth;
092
093  /** Build a Runge-Kutta integrator with the given Butcher array.
094   * @param name name of the method
095   * @param fsal indicate that the method is an <i>fsal</i>
096   * @param c time steps from Butcher array (without the first zero)
097   * @param a internal weights from Butcher array (without the first empty row)
098   * @param b propagation weights for the high order method from Butcher array
099   * @param prototype prototype of the step interpolator to use
100   * @param minStep minimal step (sign is irrelevant, regardless of
101   * integration direction, forward or backward), the last step can
102   * be smaller than this
103   * @param maxStep maximal step (sign is irrelevant, regardless of
104   * integration direction, forward or backward), the last step can
105   * be smaller than this
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 void integrate(final ExpandableStatefulODE equations, final double t)
193      throws NumberIsTooSmallException, DimensionMismatchException,
194             MaxCountExceededException, NoBracketingException {
195
196    sanityChecks(equations, t);
197    setEquations(equations);
198    final boolean forward = t > equations.getTime();
199
200    // create some internal working arrays
201    final double[] y0  = equations.getCompleteState();
202    final double[] y = y0.clone();
203    final int stages = c.length + 1;
204    final double[][] yDotK = new double[stages][y.length];
205    final double[] yTmp    = y0.clone();
206    final double[] yDotTmp = new double[y.length];
207
208    // set up an interpolator sharing the integrator arrays
209    final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
210    interpolator.reinitialize(this, yTmp, yDotK, forward,
211                              equations.getPrimaryMapper(), equations.getSecondaryMappers());
212    interpolator.storeTime(equations.getTime());
213
214    // set up integration control objects
215    stepStart         = equations.getTime();
216    double  hNew      = 0;
217    boolean firstTime = true;
218    initIntegration(equations.getTime(), y0, t);
219
220    // main integration loop
221    isLastStep = false;
222    do {
223
224      interpolator.shift();
225
226      // iterate over step size, ensuring local normalized error is smaller than 1
227      double error = 10;
228      while (error >= 1.0) {
229
230        if (firstTime || !fsal) {
231          // first stage
232          computeDerivatives(stepStart, y, yDotK[0]);
233        }
234
235        if (firstTime) {
236          final double[] scale = new double[mainSetDimension];
237          if (vecAbsoluteTolerance == null) {
238              for (int i = 0; i < scale.length; ++i) {
239                scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
240              }
241          } else {
242              for (int i = 0; i < scale.length; ++i) {
243                scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
244              }
245          }
246          hNew = initializeStep(forward, getOrder(), scale,
247                                stepStart, y, yDotK[0], yTmp, yDotK[1]);
248          firstTime = false;
249        }
250
251        stepSize = hNew;
252        if (forward) {
253            if (stepStart + stepSize >= t) {
254                stepSize = t - stepStart;
255            }
256        } else {
257            if (stepStart + stepSize <= t) {
258                stepSize = t - stepStart;
259            }
260        }
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          // reject the step and attempt to reduce error by stepsize control
290          final double factor =
291              FastMath.min(maxGrowth,
292                           FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
293          hNew = filterStep(stepSize * factor, forward, false);
294        }
295
296      }
297
298      // local error is small enough: accept the step, trigger events and step handlers
299      interpolator.storeTime(stepStart + stepSize);
300      System.arraycopy(yTmp, 0, y, 0, y0.length);
301      System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
302      stepStart = acceptStep(interpolator, y, yDotTmp, t);
303      System.arraycopy(y, 0, yTmp, 0, y.length);
304
305      if (!isLastStep) {
306
307          // prepare next step
308          interpolator.storeTime(stepStart);
309
310          if (fsal) {
311              // save the last evaluation for the next step
312              System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
313          }
314
315          // stepsize control for next step
316          final double factor =
317              FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
318          final double  scaledH    = stepSize * factor;
319          final double  nextT      = stepStart + scaledH;
320          final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
321          hNew = filterStep(scaledH, forward, nextIsLast);
322
323          final double  filteredNextT      = stepStart + hNew;
324          final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
325          if (filteredNextIsLast) {
326              hNew = t - stepStart;
327          }
328
329      }
330
331    } while (!isLastStep);
332
333    // dispatch results
334    equations.setTime(stepStart);
335    equations.setCompleteState(y);
336
337    resetInternalState();
338
339  }
340
341  /** Get the minimal reduction factor for stepsize control.
342   * @return minimal reduction factor
343   */
344  public double getMinReduction() {
345    return minReduction;
346  }
347
348  /** Set the minimal reduction factor for stepsize control.
349   * @param minReduction minimal reduction factor
350   */
351  public void setMinReduction(final double minReduction) {
352    this.minReduction = minReduction;
353  }
354
355  /** Get the maximal growth factor for stepsize control.
356   * @return maximal growth factor
357   */
358  public double getMaxGrowth() {
359    return maxGrowth;
360  }
361
362  /** Set the maximal growth factor for stepsize control.
363   * @param maxGrowth maximal growth factor
364   */
365  public void setMaxGrowth(final double maxGrowth) {
366    this.maxGrowth = maxGrowth;
367  }
368
369  /** Compute the error ratio.
370   * @param yDotK derivatives computed during the first stages
371   * @param y0 estimate of the step at the start of the step
372   * @param y1 estimate of the step at the end of the step
373   * @param h  current step
374   * @return error ratio, greater than 1 if step should be rejected
375   */
376  protected abstract double estimateError(double[][] yDotK,
377                                          double[] y0, double[] y1,
378                                          double h);
379
380}