    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
    021    import org.apache.commons.math3.exception.DimensionMismatchException;
    022    import org.apache.commons.math3.exception.MaxCountExceededException;
    023    import org.apache.commons.math3.exception.NoBracketingException;
    024    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    025    import org.apache.commons.math3.ode.AbstractIntegrator;
    026    import org.apache.commons.math3.ode.ExpandableStatefulODE;
    027    import org.apache.commons.math3.util.FastMath;
    028
    029    /**
    030     * This class implements the common part of all fixed step Runge-Kutta
    031     * integrators for Ordinary Differential Equations.
    032     *
    033     * <p>These methods are explicit Runge-Kutta methods, 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     * </pre>
    044     * </p>
    045     *
    046     * @see EulerIntegrator
    047     * @see ClassicalRungeKuttaIntegrator
    048     * @see GillIntegrator
    049     * @see MidpointIntegrator
    050     * @version $Id: RungeKuttaIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
    051     * @since 1.2
    052     */
    053
    054    public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
    055
    056        /** Time steps from Butcher array (without the first zero). */
    057        private final double[] c;
    058
    059        /** Internal weights from Butcher array (without the first empty row). */
    060        private final double[][] a;
    061
    062        /** External weights for the high order method from Butcher array. */
    063        private final double[] b;
    064
    065        /** Prototype of the step interpolator. */
    066        private final RungeKuttaStepInterpolator prototype;
    067
    068        /** Integration step. */
    069        private final double step;
    070
    071      /** Simple constructor.
    072       * Build a Runge-Kutta integrator with the given
    073       * step. The default step handler does nothing.
    074       * @param name name of the method
    075       * @param c time steps from Butcher array (without the first zero)
    076       * @param a internal weights from Butcher array (without the first empty row)
    077       * @param b propagation weights for the high order method from Butcher array
    078       * @param prototype prototype of the step interpolator to use
    079       * @param step integration step
    080       */
    081      protected RungeKuttaIntegrator(final String name,
    082                                     final double[] c, final double[][] a, final double[] b,
    083                                     final RungeKuttaStepInterpolator prototype,
    084                                     final double step) {
    085        super(name);
    086        this.c          = c;
    087        this.a          = a;
    088        this.b          = b;
    089        this.prototype  = prototype;
    090        this.step       = FastMath.abs(step);
    091      }
    092
    093      /** {@inheritDoc} */
    094      @Override
    095      public void integrate(final ExpandableStatefulODE equations, final double t)
    096          throws NumberIsTooSmallException, DimensionMismatchException,
    097                 MaxCountExceededException, NoBracketingException {
    098
    099        sanityChecks(equations, t);
    100        setEquations(equations);
    101        final boolean forward = t > equations.getTime();
    102
    103        // create some internal working arrays
    104        final double[] y0      = equations.getCompleteState();
    105        final double[] y       = y0.clone();
    106        final int stages       = c.length + 1;
    107        final double[][] yDotK = new double[stages][];
    108        for (int i = 0; i < stages; ++i) {
    109          yDotK [i] = new double[y0.length];
    110        }
    111        final double[] yTmp    = y0.clone();
    112        final double[] yDotTmp = new double[y0.length];
    113
    114        // set up an interpolator sharing the integrator arrays
    115        final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
    116        interpolator.reinitialize(this, yTmp, yDotK, forward,
    117                                  equations.getPrimaryMapper(), equations.getSecondaryMappers());
    118        interpolator.storeTime(equations.getTime());
    119
    120        // set up integration control objects
    121        stepStart = equations.getTime();
    122        stepSize  = forward ? step : -step;
    123        initIntegration(equations.getTime(), y0, t);
    124
    125        // main integration loop
    126        isLastStep = false;
    127        do {
    128
    129          interpolator.shift();
    130
    131          // first stage
    132          computeDerivatives(stepStart, y, yDotK[0]);
    133
    134          // next stages
    135          for (int k = 1; k < stages; ++k) {
    136
    137              for (int j = 0; j < y0.length; ++j) {
    138                  double sum = a[k-1][0] * yDotK[0][j];
    139                  for (int l = 1; l < k; ++l) {
    140                      sum += a[k-1][l] * yDotK[l][j];
    141                  }
    142                  yTmp[j] = y[j] + stepSize * sum;
    143              }
    144
    145              computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
    146
    147          }
    148
    149          // estimate the state at the end of the step
    150          for (int j = 0; j < y0.length; ++j) {
    151              double sum    = b[0] * yDotK[0][j];
    152              for (int l = 1; l < stages; ++l) {
    153                  sum    += b[l] * yDotK[l][j];
    154              }
    155              yTmp[j] = y[j] + stepSize * sum;
    156          }
    157
    158          // discrete events handling
    159          interpolator.storeTime(stepStart + stepSize);
    160          System.arraycopy(yTmp, 0, y, 0, y0.length);
    161          System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
    162          stepStart = acceptStep(interpolator, y, yDotTmp, t);
    163
    164          if (!isLastStep) {
    165
    166              // prepare next step
    167              interpolator.storeTime(stepStart);
    168
    169              // stepsize control for next step
    170              final double  nextT      = stepStart + stepSize;
    171              final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
    172              if (nextIsLast) {
    173                  stepSize = t - stepStart;
    174              }
    175          }
    176
    177        } while (!isLastStep);
    178
    179        // dispatch results
    180        equations.setTime(stepStart);
    181        equations.setCompleteState(y);
    182
    183        stepStart = Double.NaN;
    184        stepSize  = Double.NaN;
    185
    186      }
    187
    188    }