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: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $
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       = Math.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    
116        // set up an interpolator sharing the integrator arrays
117        AbstractStepInterpolator interpolator;
118        if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
119          final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
120          rki.reinitialize(this, yTmp, yDotK, forward);
121          interpolator = rki;
122        } else {
123          interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
124        }
125        interpolator.storeTime(t0);
126    
127        // set up integration control objects
128        stepStart = t0;
129        stepSize  = forward ? step : -step;
130        for (StepHandler handler : stepHandlers) {
131            handler.reset();
132        }
133        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
134        boolean lastStep = false;
135    
136        // main integration loop
137        while (!lastStep) {
138    
139          interpolator.shift();
140    
141          for (boolean loop = true; loop;) {
142    
143            // first stage
144            computeDerivatives(stepStart, y, yDotK[0]);
145    
146            // next stages
147            for (int k = 1; k < stages; ++k) {
148    
149              for (int j = 0; j < y0.length; ++j) {
150                double sum = a[k-1][0] * yDotK[0][j];
151                for (int l = 1; l < k; ++l) {
152                  sum += a[k-1][l] * yDotK[l][j];
153                }
154                yTmp[j] = y[j] + stepSize * sum;
155              }
156    
157              computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
158    
159            }
160    
161            // estimate the state at the end of the step
162            for (int j = 0; j < y0.length; ++j) {
163              double sum    = b[0] * yDotK[0][j];
164              for (int l = 1; l < stages; ++l) {
165                sum    += b[l] * yDotK[l][j];
166              }
167              yTmp[j] = y[j] + stepSize * sum;
168            }
169    
170            // discrete events handling
171            interpolator.storeTime(stepStart + stepSize);
172            if (manager.evaluateStep(interpolator)) {
173                final double dt = manager.getEventTime() - stepStart;
174                if (Math.abs(dt) <= Math.ulp(stepStart)) {
175                    // we cannot simply truncate the step, reject the current computation
176                    // and let the loop compute another state with the truncated step.
177                    // it is so small (much probably exactly 0 due to limited accuracy)
178                    // that the code above would fail handling it.
179                    // So we set up an artificial 0 size step by copying states
180                    interpolator.storeTime(stepStart);
181                    System.arraycopy(y, 0, yTmp, 0, y0.length);
182                    stepSize = 0;
183                    loop     = false;
184                } else {
185                    // reject the step to match exactly the next switch time
186                    stepSize = dt;
187                }
188            } else {
189              loop = false;
190            }
191    
192          }
193    
194          // the step has been accepted
195          final double nextStep = stepStart + stepSize;
196          System.arraycopy(yTmp, 0, y, 0, y0.length);
197          manager.stepAccepted(nextStep, y);
198          lastStep = manager.stop();
199    
200          // provide the step data to the step handler
201          interpolator.storeTime(nextStep);
202          for (StepHandler handler : stepHandlers) {
203              handler.handleStep(interpolator, lastStep);
204          }
205          stepStart = nextStep;
206    
207          if (manager.reset(stepStart, y) && ! lastStep) {
208            // some events handler has triggered changes that
209            // invalidate the derivatives, we need to recompute them
210            computeDerivatives(stepStart, y, yDotK[0]);
211          }
212    
213          // make sure step size is set to default before next step
214          stepSize = forward ? step : -step;
215    
216        }
217    
218        final double stopTime = stepStart;
219        stepStart = Double.NaN;
220        stepSize  = Double.NaN;
221        return stopTime;
222    
223      }
224    
225    }