1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package org.apache.commons.math.ode.nonstiff;
19
20
21import org.apache.commons.math.ode.AbstractIntegrator;
22import org.apache.commons.math.ode.DerivativeException;
23import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
24import org.apache.commons.math.ode.IntegratorException;
25import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
26import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
27import org.apache.commons.math.ode.sampling.StepHandler;
28import org.apache.commons.math.util.FastMath;
29
30/**
31 * This class implements the common part of all fixed step Runge-Kutta
32 * integrators for Ordinary Differential Equations.
33 *
34 * <p>These methods are explicit Runge-Kutta methods, their Butcher
35 * arrays are as follows :
36 * <pre>
37 *    0  |
38 *   c2  | a21
39 *   c3  | a31  a32
40 *   ... |        ...
41 *   cs  | as1  as2  ...  ass-1
42 *       |--------------------------
43 *       |  b1   b2  ...   bs-1  bs
44 * </pre>
45 * </p>
46 *
47 * @see EulerIntegrator
48 * @see ClassicalRungeKuttaIntegrator
49 * @see GillIntegrator
50 * @see MidpointIntegrator
51 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
52 * @since 1.2
53 */
54
55public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
56
57    /** Time steps from Butcher array (without the first zero). */
58    private final double[] c;
59
60    /** Internal weights from Butcher array (without the first empty row). */
61    private final double[][] a;
62
63    /** External weights for the high order method from Butcher array. */
64    private final double[] b;
65
66    /** Prototype of the step interpolator. */
67    private final RungeKuttaStepInterpolator prototype;
68
69    /** Integration step. */
70    private final double step;
71
72  /** Simple constructor.
73   * Build a Runge-Kutta integrator with the given
74   * step. The default step handler does nothing.
75   * @param name name of the method
76   * @param c time steps from Butcher array (without the first zero)
77   * @param a internal weights from Butcher array (without the first empty row)
78   * @param b propagation weights for the high order method from Butcher array
79   * @param prototype prototype of the step interpolator to use
80   * @param step integration step
81   */
82  protected RungeKuttaIntegrator(final String name,
83                                 final double[] c, final double[][] a, final double[] b,
84                                 final RungeKuttaStepInterpolator prototype,
85                                 final double step) {
86    super(name);
87    this.c          = c;
88    this.a          = a;
89    this.b          = b;
90    this.prototype  = prototype;
91    this.step       = FastMath.abs(step);
92  }
93
94  /** {@inheritDoc} */
95  public double integrate(final FirstOrderDifferentialEquations equations,
96                          final double t0, final double[] y0,
97                          final double t, final double[] y)
98  throws DerivativeException, IntegratorException {
99
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}
198