    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 org.apache.commons.math.ode.DerivativeException;
    021    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
    022    import org.apache.commons.math.ode.IntegratorException;
    023    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
    024    import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
    025    import org.apache.commons.math.ode.sampling.StepHandler;
    026    import org.apache.commons.math.util.FastMath;
    027
    028    /**
    029     * This class implements the common part of all embedded Runge-Kutta
    030     * integrators for Ordinary Differential Equations.
    031     *
    032     * <p>These methods are embedded explicit Runge-Kutta methods with two
    033     * sets of coefficients allowing to estimate the error, their Butcher
    034     * arrays are as follows :
    035     * <pre>
    036     *    0  |
    037     *   c2  | a21
    038     *   c3  | a31  a32
    039     *   ... |        ...
    040     *   cs  | as1  as2  ...  ass-1
    041     *       |--------------------------
    042     *       |  b1   b2  ...   bs-1  bs
    043     *       |  b'1  b'2 ...   b's-1 b's
    044     * </pre>
    045     * </p>
    046     *
    047     * <p>In fact, we rather use the array defined by ej = bj - b'j to
    048     * compute directly the error rather than computing two estimates and
    049     * then comparing them.</p>
    050     *
    051     * <p>Some methods are qualified as <i>fsal</i> (first same as last)
    052     * methods. This means the last evaluation of the derivatives in one
    053     * step is the same as the first in the next step. Then, this
    054     * evaluation can be reused from one step to the next one and the cost
    055     * of such a method is really s-1 evaluations despite the method still
    056     * has s stages. This behaviour is true only for successful steps, if
    057     * the step is rejected after the error estimation phase, no
    058     * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
    059     * asi = bi for all i.</p>
    060     *
    061     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    062     * @since 1.2
    063     */
    064
    065    public abstract class EmbeddedRungeKuttaIntegrator
    066      extends AdaptiveStepsizeIntegrator {
    067
    068        /** Indicator for <i>fsal</i> methods. */
    069        private final boolean fsal;
    070
    071        /** Time steps from Butcher array (without the first zero). */
    072        private final double[] c;
    073
    074        /** Internal weights from Butcher array (without the first empty row). */
    075        private final double[][] a;
    076
    077        /** External weights for the high order method from Butcher array. */
    078        private final double[] b;
    079
    080        /** Prototype of the step interpolator. */
    081        private final RungeKuttaStepInterpolator prototype;
    082
    083        /** Stepsize control exponent. */
    084        private final double exp;
    085
    086        /** Safety factor for stepsize control. */
    087        private double safety;
    088
    089        /** Minimal reduction factor for stepsize control. */
    090        private double minReduction;
    091
    092        /** Maximal growth factor for stepsize control. */
    093        private double maxGrowth;
    094
    095      /** Build a Runge-Kutta integrator with the given Butcher array.
    096       * @param name name of the method
    097       * @param fsal indicate that the method is an <i>fsal</i>
    098       * @param c time steps from Butcher array (without the first zero)
    099       * @param a internal weights from Butcher array (without the first empty row)
    100       * @param b propagation weights for the high order method from Butcher array
    101       * @param prototype prototype of the step interpolator to use
    102       * @param minStep minimal step (must be positive even for backward
    103       * integration), the last step can be smaller than this
    104       * @param maxStep maximal step (must be positive even for backward
    105       * integration)
    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 double integrate(final FirstOrderDifferentialEquations equations,
    193                              final double t0, final double[] y0,
    194                              final double t, final double[] y)
    195      throws DerivativeException, IntegratorException {
    196
    197        sanityChecks(equations, t0, y0, t, y);
    198        setEquations(equations);
    199        resetEvaluations();
    200        final boolean forward = t > t0;
    201
    202        // create some internal working arrays
    203        final int stages = c.length + 1;
    204        if (y != y0) {
    205          System.arraycopy(y0, 0, y, 0, y0.length);
    206        }
    207        final double[][] yDotK = new double[stages][y0.length];
    208        final double[] yTmp    = new double[y0.length];
    209        final double[] yDotTmp = new double[y0.length];
    210
    211        // set up an interpolator sharing the integrator arrays
    212        AbstractStepInterpolator interpolator;
    213        if (requiresDenseOutput()) {
    214          final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
    215          rki.reinitialize(this, yTmp, yDotK, forward);
    216          interpolator = rki;
    217        } else {
    218          interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
    219        }
    220        interpolator.storeTime(t0);
    221
    222        // set up integration control objects
    223        stepStart         = t0;
    224        double  hNew      = 0;
    225        boolean firstTime = true;
    226        for (StepHandler handler : stepHandlers) {
    227            handler.reset();
    228        }
    229        setStateInitialized(false);
    230
    231        // main integration loop
    232        isLastStep = false;
    233        do {
    234
    235          interpolator.shift();
    236
    237          // iterate over step size, ensuring local normalized error is smaller than 1
    238          double error = 10;
    239          while (error >= 1.0) {
    240
    241            if (firstTime || !fsal) {
    242              // first stage
    243              computeDerivatives(stepStart, y, yDotK[0]);
    244            }
    245
    246            if (firstTime) {
    247              final double[] scale = new double[mainSetDimension];
    248              if (vecAbsoluteTolerance == null) {
    249                  for (int i = 0; i < scale.length; ++i) {
    250                    scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
    251                  }
    252              } else {
    253                  for (int i = 0; i < scale.length; ++i) {
    254                    scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
    255                  }
    256              }
    257              hNew = initializeStep(equations, forward, getOrder(), scale,
    258                                    stepStart, y, yDotK[0], yTmp, yDotK[1]);
    259              firstTime = false;
    260            }
    261
    262            stepSize = hNew;
    263
    264            // next stages
    265            for (int k = 1; k < stages; ++k) {
    266
    267              for (int j = 0; j < y0.length; ++j) {
    268                double sum = a[k-1][0] * yDotK[0][j];
    269                for (int l = 1; l < k; ++l) {
    270                  sum += a[k-1][l] * yDotK[l][j];
    271                }
    272                yTmp[j] = y[j] + stepSize * sum;
    273              }
    274
    275              computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
    276
    277            }
    278
    279            // estimate the state at the end of the step
    280            for (int j = 0; j < y0.length; ++j) {
    281              double sum    = b[0] * yDotK[0][j];
    282              for (int l = 1; l < stages; ++l) {
    283                sum    += b[l] * yDotK[l][j];
    284              }
    285              yTmp[j] = y[j] + stepSize * sum;
    286            }
    287
    288            // estimate the error at the end of the step
    289            error = estimateError(yDotK, y, yTmp, stepSize);
    290            if (error >= 1.0) {
    291              // reject the step and attempt to reduce error by stepsize control
    292              final double factor =
    293                  FastMath.min(maxGrowth,
    294                               FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
    295              hNew = filterStep(stepSize * factor, forward, false);
    296            }
    297
    298          }
    299
    300          // local error is small enough: accept the step, trigger events and step handlers
    301          interpolator.storeTime(stepStart + stepSize);
    302          System.arraycopy(yTmp, 0, y, 0, y0.length);
    303          System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
    304          stepStart = acceptStep(interpolator, y, yDotTmp, t);
    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        final double stopTime = stepStart;
    335        resetInternalState();
    336        return stopTime;
    337
    338      }
    339
    340      /** Get the minimal reduction factor for stepsize control.
    341       * @return minimal reduction factor
    342       */
    343      public double getMinReduction() {
    344        return minReduction;
    345      }
    346
    347      /** Set the minimal reduction factor for stepsize control.
    348       * @param minReduction minimal reduction factor
    349       */
    350      public void setMinReduction(final double minReduction) {
    351        this.minReduction = minReduction;
    352      }
    353
    354      /** Get the maximal growth factor for stepsize control.
    355       * @return maximal growth factor
    356       */
    357      public double getMaxGrowth() {
    358        return maxGrowth;
    359      }
    360
    361      /** Set the maximal growth factor for stepsize control.
    362       * @param maxGrowth maximal growth factor
    363       */
    364      public void setMaxGrowth(final double maxGrowth) {
    365        this.maxGrowth = maxGrowth;
    366      }
    367
    368      /** Compute the error ratio.
    369       * @param yDotK derivatives computed during the first stages
    370       * @param y0 estimate of the step at the start of the step
    371       * @param y1 estimate of the step at the end of the step
    372       * @param h  current step
    373       * @return error ratio, greater than 1 if step should be rejected
    374       */
    375      protected abstract double estimateError(double[][] yDotK,
    376                                              double[] y0, double[] y1,
    377                                              double h);
    378
    379    }