/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.ode.AbstractIntegrator;
import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.util.Cloner;
import org.apache.commons.math3.util.FastMath;
/**
* This class implements the common part of all fixed step Runge-Kutta
* integrators for Ordinary Differential Equations.
*
* <p>These methods are explicit Runge-Kutta methods, their Butcher
* arrays are as follows :
* <pre>
* 0 |
* c2 | a21
* c3 | a31 a32
* ... | ...
* cs | as1 as2 ... ass-1
* |--------------------------
* | b1 b2 ... bs-1 bs
* </pre>
* </p>
*
* @see EulerIntegrator
* @see ClassicalRungeKuttaIntegrator
* @see GillIntegrator
* @see MidpointIntegrator
* @since 1.2
*/
public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
/** Time steps from Butcher array (without the first zero). */
private final double[] c;
/** Internal weights from Butcher array (without the first empty row). */
private final double[][] a;
/** External weights for the high order method from Butcher array. */
private final double[] b;
/** Prototype of the step interpolator. */
private final RungeKuttaStepInterpolator prototype;
/** Integration step. */
private final double step;
/** Simple constructor.
* Build a Runge-Kutta integrator with the given
* step. The default step handler does nothing.
* @param name name of the method
* @param c time steps from Butcher array (without the first zero)
* @param a internal weights from Butcher array (without the first empty row)
* @param b propagation weights for the high order method from Butcher array
* @param prototype prototype of the step interpolator to use
* @param step integration step
*/
protected RungeKuttaIntegrator(final String name,
final double[] c, final double[][] a, final double[] b,
final RungeKuttaStepInterpolator prototype,
final double step) {
super(name);
this.c = c;
this.a = a;
this.b = b;
this.prototype = prototype;
this.step = Math.abs(step);
}
/** {@inheritDoc} */
@Override
public void integrate(final ExpandableStatefulODE equations, final double t)
throws NumberIsTooSmallException, DimensionMismatchException,
MaxCountExceededException, NoBracketingException {
sanityChecks(equations, t);
setEquations(equations);
final boolean forward = t > equations.getTime();
// create some internal working arrays
final double[] y0 = equations.getCompleteState();
final double[] y = Cloner.clone(y0);
final int stages = c.length + 1;
final double[][] yDotK = new double[stages][];
for (int i = 0; i < stages; ++i) {
yDotK [i] = new double[y0.length];
}
final double[] yTmp = Cloner.clone(y0);
final double[] yDotTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
interpolator.reinitialize(this, yTmp, yDotK, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
interpolator.storeTime(equations.getTime());
// set up integration control objects
stepStart = equations.getTime();
if (forward) {
if (stepStart + step >= t) {
stepSize = t - stepStart;
} else {
stepSize = step;
}
} else {
if (stepStart - step <= t) {
stepSize = t - stepStart;
} else {
stepSize = -step;
}
}
initIntegration(equations.getTime(), y0, t);
// main integration loop
isLastStep = false;
do {
interpolator.shift();
// first stage
computeDerivatives(stepStart, y, yDotK[0]);
// next stages
for (int k = 1; k < stages; ++k) {
for (int j = 0; j < y0.length; ++j) {
double sum = a[k-1][0] * yDotK[0][j];
for (int l = 1; l < k; ++l) {
sum += a[k-1][l] * yDotK[l][j];
}
yTmp[j] = y[j] + stepSize * sum;
}
computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
}
// estimate the state at the end of the step
for (int j = 0; j < y0.length; ++j) {
double sum = b[0] * yDotK[0][j];
for (int l = 1; l < stages; ++l) {
sum += b[l] * yDotK[l][j];
}
yTmp[j] = y[j] + stepSize * sum;
}
// discrete events handling
interpolator.storeTime(stepStart + stepSize);
System.arraycopy(yTmp, 0, y, 0, y0.length);
System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
stepStart = acceptStep(interpolator, y, yDotTmp, t);
if (!isLastStep) {
// prepare next step
interpolator.storeTime(stepStart);
// stepsize control for next step
final double nextT = stepStart + stepSize;
final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
if (nextIsLast) {
stepSize = t - stepStart;
}
}
} while (!isLastStep);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
stepStart = Double.NaN;
stepSize = Double.NaN;
}
/** Fast computation of a single step of ODE integration.
* <p>This method is intended for the limited use case of
* very fast computation of only one step without using any of the
* rich features of general integrators that may take some time
* to set up (i.e. no step handlers, no events handlers, no additional
* states, no interpolators, no error control, no evaluations count,
* no sanity checks ...). It handles the strict minimum of computation,
* so it can be embedded in outer loops.</p>
* <p>
* This method is <em>not</em> used at all by the {@link #integrate(ExpandableStatefulODE, double)}
* method. It also completely ignores the step set at construction time, and
* uses only a single step to go from {@code t0} to {@code t}.
* </p>
* <p>
* As this method does not use any of the state-dependent features of the integrator,
* it should be reasonably thread-safe <em>if and only if</em> the provided differential
* equations are themselves thread-safe.
* </p>
* @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param t target time for the integration
* (can be set to a value smaller than {@code t0} for backward integration)
* @return state vector at {@code t}
*/
public double[] singleStep(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0, final double t) {
// create some internal working arrays
final double[] y = Cloner.clone(y0);
final int stages = c.length + 1;
final double[][] yDotK = new double[stages][];
for (int i = 0; i < stages; ++i) {
yDotK [i] = new double[y0.length];
}
final double[] yTmp = Cloner.clone(y0);
// first stage
final double h = t - t0;
equations.computeDerivatives(t0, y, yDotK[0]);
// next stages
for (int k = 1; k < stages; ++k) {
for (int j = 0; j < y0.length; ++j) {
double sum = a[k-1][0] * yDotK[0][j];
for (int l = 1; l < k; ++l) {
sum += a[k-1][l] * yDotK[l][j];
}
yTmp[j] = y[j] + h * sum;
}
equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]);
}
// estimate the state at the end of the step
for (int j = 0; j < y0.length; ++j) {
double sum = b[0] * yDotK[0][j];
for (int l = 1; l < stages; ++l) {
sum += b[l] * yDotK[l][j];
}
y[j] += h * sum;
}
return y;
}
}