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.legacy.ode.nonstiff;
019
020
021import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
022import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
023import org.apache.commons.math4.legacy.exception.NoBracketingException;
024import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
025import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
026import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
027import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
028import org.apache.commons.math4.core.jdkmath.JdkMath;
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       = JdkMath.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      // estimate the state at the end of the step
160      for (int j = 0; j < y0.length; ++j) {
161          double sum    = b[0] * yDotK[0][j];
162          for (int l = 1; l < stages; ++l) {
163              sum    += b[l] * yDotK[l][j];
164          }
165          yTmp[j] = y[j] + stepSize * sum;
166      }
167
168      // discrete events handling
169      interpolator.storeTime(stepStart + stepSize);
170      System.arraycopy(yTmp, 0, y, 0, y0.length);
171      System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
172      stepStart = acceptStep(interpolator, y, yDotTmp, t);
173
174      if (!isLastStep) {
175
176          // prepare next step
177          interpolator.storeTime(stepStart);
178
179          // stepsize control for next step
180          final double  nextT      = stepStart + stepSize;
181          final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
182          if (nextIsLast) {
183              stepSize = t - stepStart;
184          }
185      }
186    } while (!isLastStep);
187
188    // dispatch results
189    equations.setTime(stepStart);
190    equations.setCompleteState(y);
191
192    stepStart = Double.NaN;
193    stepSize  = Double.NaN;
194  }
195
196  /** Fast computation of a single step of ODE integration.
197   * <p>This method is intended for the limited use case of
198   * very fast computation of only one step without using any of the
199   * rich features of general integrators that may take some time
200   * to set up (i.e. no step handlers, no events handlers, no additional
201   * states, no interpolators, no error control, no evaluations count,
202   * no sanity checks ...). It handles the strict minimum of computation,
203   * so it can be embedded in outer loops.</p>
204   * <p>
205   * This method is <em>not</em> used at all by the {@link #integrate(ExpandableStatefulODE, double)}
206   * method. It also completely ignores the step set at construction time, and
207   * uses only a single step to go from {@code t0} to {@code t}.
208   * </p>
209   * <p>
210   * As this method does not use any of the state-dependent features of the integrator,
211   * it should be reasonably thread-safe <em>if and only if</em> the provided differential
212   * equations are themselves thread-safe.
213   * </p>
214   * @param equations differential equations to integrate
215   * @param t0 initial time
216   * @param y0 initial value of the state vector at t0
217   * @param t target time for the integration
218   * (can be set to a value smaller than {@code t0} for backward integration)
219   * @return state vector at {@code t}
220   */
221  public double[] singleStep(final FirstOrderDifferentialEquations equations,
222                             final double t0, final double[] y0, final double t) {
223
224      // create some internal working arrays
225      final double[] y       = y0.clone();
226      final int stages       = c.length + 1;
227      final double[][] yDotK = new double[stages][];
228      for (int i = 0; i < stages; ++i) {
229          yDotK [i] = new double[y0.length];
230      }
231      final double[] yTmp    = y0.clone();
232
233      // first stage
234      final double h = t - t0;
235      equations.computeDerivatives(t0, y, yDotK[0]);
236
237      // next stages
238      for (int k = 1; k < stages; ++k) {
239
240          for (int j = 0; j < y0.length; ++j) {
241              double sum = a[k-1][0] * yDotK[0][j];
242              for (int l = 1; l < k; ++l) {
243                  sum += a[k-1][l] * yDotK[l][j];
244              }
245              yTmp[j] = y[j] + h * sum;
246          }
247
248          equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]);
249      }
250
251      // estimate the state at the end of the step
252      for (int j = 0; j < y0.length; ++j) {
253          double sum = b[0] * yDotK[0][j];
254          for (int l = 1; l < stages; ++l) {
255              sum += b[l] * yDotK[l][j];
256          }
257          y[j] += h * sum;
258      }
259
260      return y;
261  }
262}