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
018package org.apache.commons.math4.ode.nonstiff;
019
020
021import org.apache.commons.math4.exception.DimensionMismatchException;
022import org.apache.commons.math4.exception.MaxCountExceededException;
023import org.apache.commons.math4.exception.NoBracketingException;
024import org.apache.commons.math4.exception.NumberIsTooSmallException;
025import org.apache.commons.math4.ode.AbstractIntegrator;
026import org.apache.commons.math4.ode.ExpandableStatefulODE;
027import org.apache.commons.math4.ode.FirstOrderDifferentialEquations;
028import org.apache.commons.math4.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 *
046 * @see EulerIntegrator
047 * @see ClassicalRungeKuttaIntegrator
048 * @see GillIntegrator
049 * @see MidpointIntegrator
050 * @since 1.2
051 */
052
053public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
054
055    /** Time steps from Butcher array (without the first zero). */
056    private final double[] c;
057
058    /** Internal weights from Butcher array (without the first empty row). */
059    private final double[][] a;
060
061    /** External weights for the high order method from Butcher array. */
062    private final double[] b;
063
064    /** Prototype of the step interpolator. */
065    private final RungeKuttaStepInterpolator prototype;
066
067    /** Integration step. */
068    private final double step;
069
070  /** Simple constructor.
071   * Build a Runge-Kutta integrator with the given
072   * step. The default step handler does nothing.
073   * @param name name of the method
074   * @param c time steps from Butcher array (without the first zero)
075   * @param a internal weights from Butcher array (without the first empty row)
076   * @param b propagation weights for the high order method from Butcher array
077   * @param prototype prototype of the step interpolator to use
078   * @param step integration step
079   */
080  protected RungeKuttaIntegrator(final String name,
081                                 final double[] c, final double[][] a, final double[] b,
082                                 final RungeKuttaStepInterpolator prototype,
083                                 final double step) {
084    super(name);
085    this.c          = c;
086    this.a          = a;
087    this.b          = b;
088    this.prototype  = prototype;
089    this.step       = FastMath.abs(step);
090  }
091
092  /** {@inheritDoc} */
093  @Override
094  public void integrate(final ExpandableStatefulODE equations, final double t)
095      throws NumberIsTooSmallException, DimensionMismatchException,
096             MaxCountExceededException, NoBracketingException {
097
098    sanityChecks(equations, t);
099    setEquations(equations);
100    final boolean forward = t > equations.getTime();
101
102    // create some internal working arrays
103    final double[] y0      = equations.getCompleteState();
104    final double[] y       = y0.clone();
105    final int stages       = c.length + 1;
106    final double[][] yDotK = new double[stages][];
107    for (int i = 0; i < stages; ++i) {
108      yDotK [i] = new double[y0.length];
109    }
110    final double[] yTmp    = y0.clone();
111    final double[] yDotTmp = new double[y0.length];
112
113    // set up an interpolator sharing the integrator arrays
114    final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
115    interpolator.reinitialize(this, yTmp, yDotK, forward,
116                              equations.getPrimaryMapper(), equations.getSecondaryMappers());
117    interpolator.storeTime(equations.getTime());
118
119    // set up integration control objects
120    stepStart = equations.getTime();
121    if (forward) {
122        if (stepStart + step >= t) {
123            stepSize = t - stepStart;
124        } else {
125            stepSize = step;
126        }
127    } else {
128        if (stepStart - step <= t) {
129            stepSize = t - stepStart;
130        } else {
131            stepSize = -step;
132        }
133    }
134    initIntegration(equations.getTime(), y0, t);
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    // dispatch results
191    equations.setTime(stepStart);
192    equations.setCompleteState(y);
193
194    stepStart = Double.NaN;
195    stepSize  = Double.NaN;
196
197  }
198
199  /** Fast computation of a single step of ODE integration.
200   * <p>This method is intended for the limited use case of
201   * very fast computation of only one step without using any of the
202   * rich features of general integrators that may take some time
203   * to set up (i.e. no step handlers, no events handlers, no additional
204   * states, no interpolators, no error control, no evaluations count,
205   * no sanity checks ...). It handles the strict minimum of computation,
206   * so it can be embedded in outer loops.</p>
207   * <p>
208   * This method is <em>not</em> used at all by the {@link #integrate(ExpandableStatefulODE, double)}
209   * method. It also completely ignores the step set at construction time, and
210   * uses only a single step to go from {@code t0} to {@code t}.
211   * </p>
212   * <p>
213   * As this method does not use any of the state-dependent features of the integrator,
214   * it should be reasonably thread-safe <em>if and only if</em> the provided differential
215   * equations are themselves thread-safe.
216   * </p>
217   * @param equations differential equations to integrate
218   * @param t0 initial time
219   * @param y0 initial value of the state vector at t0
220   * @param t target time for the integration
221   * (can be set to a value smaller than {@code t0} for backward integration)
222   * @return state vector at {@code t}
223   */
224  public double[] singleStep(final FirstOrderDifferentialEquations equations,
225                             final double t0, final double[] y0, final double t) {
226
227      // create some internal working arrays
228      final double[] y       = y0.clone();
229      final int stages       = c.length + 1;
230      final double[][] yDotK = new double[stages][];
231      for (int i = 0; i < stages; ++i) {
232          yDotK [i] = new double[y0.length];
233      }
234      final double[] yTmp    = y0.clone();
235
236      // first stage
237      final double h = t - t0;
238      equations.computeDerivatives(t0, y, yDotK[0]);
239
240      // next stages
241      for (int k = 1; k < stages; ++k) {
242
243          for (int j = 0; j < y0.length; ++j) {
244              double sum = a[k-1][0] * yDotK[0][j];
245              for (int l = 1; l < k; ++l) {
246                  sum += a[k-1][l] * yDotK[l][j];
247              }
248              yTmp[j] = y[j] + h * sum;
249          }
250
251          equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]);
252
253      }
254
255      // estimate the state at the end of the step
256      for (int j = 0; j < y0.length; ++j) {
257          double sum = b[0] * yDotK[0][j];
258          for (int l = 1; l < stages; ++l) {
259              sum += b[l] * yDotK[l][j];
260          }
261          y[j] += h * sum;
262      }
263
264      return y;
265
266  }
267
268}