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.events.CombinedEventsManager;
026    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
027    import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
028    import org.apache.commons.math.ode.sampling.StepHandler;
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: 785473 $ $Date: 2009-06-17 00:02:35 -0400 (Wed, 17 Jun 2009) $
052     * @since 1.2
053     */
054    
055    public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
056    
057      /** Simple constructor.
058       * Build a Runge-Kutta integrator with the given
059       * step. The default step handler does nothing.
060       * @param name name of the method
061       * @param c time steps from Butcher array (without the first zero)
062       * @param a internal weights from Butcher array (without the first empty row)
063       * @param b propagation weights for the high order method from Butcher array
064       * @param prototype prototype of the step interpolator to use
065       * @param step integration step
066       */
067      protected RungeKuttaIntegrator(final String name,
068                                     final double[] c, final double[][] a, final double[] b,
069                                     final RungeKuttaStepInterpolator prototype,
070                                     final double step) {
071        super(name);
072        this.c          = c;
073        this.a          = a;
074        this.b          = b;
075        this.prototype  = prototype;
076        this.step       = Math.abs(step);
077      }
078    
079      /** {@inheritDoc} */
080      public double integrate(final FirstOrderDifferentialEquations equations,
081                              final double t0, final double[] y0,
082                              final double t, final double[] y)
083      throws DerivativeException, IntegratorException {
084    
085        sanityChecks(equations, t0, y0, t, y);
086        setEquations(equations);
087        resetEvaluations();
088        final boolean forward = (t > t0);
089    
090        // create some internal working arrays
091        final int stages = c.length + 1;
092        if (y != y0) {
093          System.arraycopy(y0, 0, y, 0, y0.length);
094        }
095        final double[][] yDotK = new double[stages][];
096        for (int i = 0; i < stages; ++i) {
097          yDotK [i] = new double[y0.length];
098        }
099        final double[] yTmp = new double[y0.length];
100    
101        // set up an interpolator sharing the integrator arrays
102        AbstractStepInterpolator interpolator;
103        if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
104          final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
105          rki.reinitialize(this, yTmp, yDotK, forward);
106          interpolator = rki;
107        } else {
108          interpolator = new DummyStepInterpolator(yTmp, forward);
109        }
110        interpolator.storeTime(t0);
111    
112        // set up integration control objects
113        stepStart = t0;
114        stepSize  = forward ? step : -step;
115        for (StepHandler handler : stepHandlers) {
116            handler.reset();
117        }
118        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
119        boolean lastStep = false;
120    
121        // main integration loop
122        while (!lastStep) {
123    
124          interpolator.shift();
125    
126          for (boolean loop = true; loop;) {
127    
128            // first stage
129            computeDerivatives(stepStart, y, yDotK[0]);
130    
131            // next stages
132            for (int k = 1; k < stages; ++k) {
133    
134              for (int j = 0; j < y0.length; ++j) {
135                double sum = a[k-1][0] * yDotK[0][j];
136                for (int l = 1; l < k; ++l) {
137                  sum += a[k-1][l] * yDotK[l][j];
138                }
139                yTmp[j] = y[j] + stepSize * sum;
140              }
141    
142              computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
143    
144            }
145    
146            // estimate the state at the end of the step
147            for (int j = 0; j < y0.length; ++j) {
148              double sum    = b[0] * yDotK[0][j];
149              for (int l = 1; l < stages; ++l) {
150                sum    += b[l] * yDotK[l][j];
151              }
152              yTmp[j] = y[j] + stepSize * sum;
153            }
154    
155            // discrete events handling
156            interpolator.storeTime(stepStart + stepSize);
157            if (manager.evaluateStep(interpolator)) {
158                final double dt = manager.getEventTime() - stepStart;
159                if (Math.abs(dt) <= Math.ulp(stepStart)) {
160                    // rejecting the step would lead to a too small next step, we accept it
161                    loop = false;
162                } else {
163                    // reject the step to match exactly the next switch time
164                    stepSize = dt;
165                }
166            } else {
167              loop = false;
168            }
169    
170          }
171    
172          // the step has been accepted
173          final double nextStep = stepStart + stepSize;
174          System.arraycopy(yTmp, 0, y, 0, y0.length);
175          manager.stepAccepted(nextStep, y);
176          lastStep = manager.stop();
177    
178          // provide the step data to the step handler
179          interpolator.storeTime(nextStep);
180          for (StepHandler handler : stepHandlers) {
181              handler.handleStep(interpolator, lastStep);
182          }
183          stepStart = nextStep;
184    
185          if (manager.reset(stepStart, y) && ! lastStep) {
186            // some events handler has triggered changes that
187            // invalidate the derivatives, we need to recompute them
188            computeDerivatives(stepStart, y, yDotK[0]);
189          }
190    
191          // make sure step size is set to default before next step
192          stepSize = forward ? step : -step;
193    
194        }
195    
196        final double stopTime = stepStart;
197        stepStart = Double.NaN;
198        stepSize  = Double.NaN;
199        return stopTime;
200    
201      }
202    
203      /** Time steps from Butcher array (without the first zero). */
204      private double[] c;
205    
206      /** Internal weights from Butcher array (without the first empty row). */
207      private double[][] a;
208    
209      /** External weights for the high order method from Butcher array. */
210      private double[] b;
211    
212      /** Prototype of the step interpolator. */
213      private RungeKuttaStepInterpolator prototype;
214                                             
215      /** Integration step. */
216      private double step;
217    
218    }