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    }