    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
    021    import org.apache.commons.math.ode.AbstractIntegrator;
    022    import org.apache.commons.math.ode.DerivativeException;
    023    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
    024    import org.apache.commons.math.ode.IntegratorException;
    025    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
    026    import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
    027    import org.apache.commons.math.ode.sampling.StepHandler;
    028    import org.apache.commons.math.util.FastMath;
    029
    030    /**
    031     * This class implements the common part of all fixed step Runge-Kutta
    032     * integrators for Ordinary Differential Equations.
    033     *
    034     * <p>These methods are explicit Runge-Kutta methods, their Butcher
    035     * arrays are as follows :
    036     * <pre>
    037     *    0  |
    038     *   c2  | a21
    039     *   c3  | a31  a32
    040     *   ... |        ...
    041     *   cs  | as1  as2  ...  ass-1
    042     *       |--------------------------
    043     *       |  b1   b2  ...   bs-1  bs
    044     * </pre>
    045     * </p>
    046     *
    047     * @see EulerIntegrator
    048     * @see ClassicalRungeKuttaIntegrator
    049     * @see GillIntegrator
    050     * @see MidpointIntegrator
    051     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    052     * @since 1.2
    053     */
    054
    055    public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
    056
    057        /** Time steps from Butcher array (without the first zero). */
    058        private final double[] c;
    059
    060        /** Internal weights from Butcher array (without the first empty row). */
    061        private final double[][] a;
    062
    063        /** External weights for the high order method from Butcher array. */
    064        private final double[] b;
    065
    066        /** Prototype of the step interpolator. */
    067        private final RungeKuttaStepInterpolator prototype;
    068
    069        /** Integration step. */
    070        private final double step;
    071
    072      /** Simple constructor.
    073       * Build a Runge-Kutta integrator with the given
    074       * step. The default step handler does nothing.
    075       * @param name name of the method
    076       * @param c time steps from Butcher array (without the first zero)
    077       * @param a internal weights from Butcher array (without the first empty row)
    078       * @param b propagation weights for the high order method from Butcher array
    079       * @param prototype prototype of the step interpolator to use
    080       * @param step integration step
    081       */
    082      protected RungeKuttaIntegrator(final String name,
    083                                     final double[] c, final double[][] a, final double[] b,
    084                                     final RungeKuttaStepInterpolator prototype,
    085                                     final double step) {
    086        super(name);
    087        this.c          = c;
    088        this.a          = a;
    089        this.b          = b;
    090        this.prototype  = prototype;
    091        this.step       = FastMath.abs(step);
    092      }
    093
    094      /** {@inheritDoc} */
    095      public double integrate(final FirstOrderDifferentialEquations equations,
    096                              final double t0, final double[] y0,
    097                              final double t, final double[] y)
    098      throws DerivativeException, IntegratorException {
    099
    100        sanityChecks(equations, t0, y0, t, y);
    101        setEquations(equations);
    102        resetEvaluations();
    103        final boolean forward = t > t0;
    104
    105        // create some internal working arrays
    106        final int stages = c.length + 1;
    107        if (y != y0) {
    108          System.arraycopy(y0, 0, y, 0, y0.length);
    109        }
    110        final double[][] yDotK = new double[stages][];
    111        for (int i = 0; i < stages; ++i) {
    112          yDotK [i] = new double[y0.length];
    113        }
    114        final double[] yTmp    = new double[y0.length];
    115        final double[] yDotTmp = new double[y0.length];
    116
    117        // set up an interpolator sharing the integrator arrays
    118        AbstractStepInterpolator interpolator;
    119        if (requiresDenseOutput()) {
    120          final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
    121          rki.reinitialize(this, yTmp, yDotK, forward);
    122          interpolator = rki;
    123        } else {
    124          interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
    125        }
    126        interpolator.storeTime(t0);
    127
    128        // set up integration control objects
    129        stepStart = t0;
    130        stepSize  = forward ? step : -step;
    131        for (StepHandler handler : stepHandlers) {
    132            handler.reset();
    133        }
    134        setStateInitialized(false);
    135
    136        // main integration loop
    137        isLastStep = false;
    138        do {
    139
    140          interpolator.shift();
    141
    142          // first stage
    143          computeDerivatives(stepStart, y, yDotK[0]);
    144
    145          // next stages
    146          for (int k = 1; k < stages; ++k) {
    147
    148              for (int j = 0; j < y0.length; ++j) {
    149                  double sum = a[k-1][0] * yDotK[0][j];
    150                  for (int l = 1; l < k; ++l) {
    151                      sum += a[k-1][l] * yDotK[l][j];
    152                  }
    153                  yTmp[j] = y[j] + stepSize * sum;
    154              }
    155
    156              computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
    157
    158          }
    159
    160          // estimate the state at the end of the step
    161          for (int j = 0; j < y0.length; ++j) {
    162              double sum    = b[0] * yDotK[0][j];
    163              for (int l = 1; l < stages; ++l) {
    164                  sum    += b[l] * yDotK[l][j];
    165              }
    166              yTmp[j] = y[j] + stepSize * sum;
    167          }
    168
    169          // discrete events handling
    170          interpolator.storeTime(stepStart + stepSize);
    171          System.arraycopy(yTmp, 0, y, 0, y0.length);
    172          System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
    173          stepStart = acceptStep(interpolator, y, yDotTmp, t);
    174
    175          if (!isLastStep) {
    176
    177              // prepare next step
    178              interpolator.storeTime(stepStart);
    179
    180              // stepsize control for next step
    181              final double  nextT      = stepStart + stepSize;
    182              final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
    183              if (nextIsLast) {
    184                  stepSize = t - stepStart;
    185              }
    186          }
    187
    188        } while (!isLastStep);
    189
    190        final double stopTime = stepStart;
    191        stepStart = Double.NaN;
    192        stepSize  = Double.NaN;
    193        return stopTime;
    194
    195      }
    196
    197    }