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