    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    }