RungeKuttaIntegrator.java
/*
* Licensed to the Hipparchus project under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The Hipparchus project 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
*
* https://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.hipparchus.ode.nonstiff;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.ode.AbstractIntegrator;
import org.hipparchus.ode.EquationsMapper;
import org.hipparchus.ode.ExpandableODE;
import org.hipparchus.ode.LocalizedODEFormats;
import org.hipparchus.ode.ODEState;
import org.hipparchus.ode.ODEStateAndDerivative;
import org.hipparchus.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 :</p>
* <pre>
* 0 |
* c2 | a21
* c3 | a31 a32
* ... | ...
* cs | as1 as2 ... ass-1
* |--------------------------
* | b1 b2 ... bs-1 bs
* </pre>
*
* @see EulerIntegrator
* @see ClassicalRungeKuttaIntegrator
* @see GillIntegrator
* @see MidpointIntegrator
*/
public abstract class RungeKuttaIntegrator extends AbstractIntegrator implements ExplicitRungeKuttaIntegrator {
/** 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;
/** 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 step integration step
*/
protected RungeKuttaIntegrator(final String name, final double step) {
super(name);
this.c = getC();
this.a = getA();
this.b = getB();
this.step = FastMath.abs(step);
}
/** Getter for the default, positive step-size assigned at constructor level.
* @return step
*/
public double getDefaultStep() {
return this.step;
}
/** Create an interpolator.
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param mapper equations mapper for the all equations
* @return external weights for the high order method from Butcher array
*/
protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
ODEStateAndDerivative globalPreviousState,
ODEStateAndDerivative globalCurrentState,
EquationsMapper mapper);
/** {@inheritDoc} */
@Override
public ODEStateAndDerivative integrate(final ExpandableODE equations,
final ODEState initialState, final double finalTime)
throws MathIllegalArgumentException, MathIllegalStateException {
sanityChecks(initialState, finalTime);
setStepStart(initIntegration(equations, initialState, finalTime));
final boolean forward = finalTime > initialState.getTime();
// create some internal working arrays
final int stages = c.length + 1;
double[] y = getStepStart().getCompleteState();
final double[][] yDotK = new double[stages][];
final double[] yTmp = new double[y.length];
// set up integration control objects
if (forward) {
if (getStepStart().getTime() + step >= finalTime) {
setStepSize(finalTime - getStepStart().getTime());
} else {
setStepSize(step);
}
} else {
if (getStepStart().getTime() - step <= finalTime) {
setStepSize(finalTime - getStepStart().getTime());
} else {
setStepSize(-step);
}
}
// main integration loop
setIsLastStep(false);
do {
// first stage
y = getStepStart().getCompleteState();
yDotK[0] = getStepStart().getCompleteDerivative();
// next stages
ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(), y,
getStepSize(), a, c, yDotK);
incrementEvaluations(stages - 1);
// estimate the state at the end of the step
for (int j = 0; j < y.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] + getStepSize() * sum;
if (Double.isNaN(yTmp[j])) {
throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
getStepStart().getTime() + getStepSize());
}
}
final double stepEnd = getStepStart().getTime() + getStepSize();
final double[] yDotTmp = computeDerivatives(stepEnd, yTmp);
final ODEStateAndDerivative stateTmp =
equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
// discrete events handling
System.arraycopy(yTmp, 0, y, 0, y.length);
setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp,
equations.getMapper()),
finalTime));
if (!isLastStep()) {
// stepsize control for next step
final double nextT = getStepStart().getTime() + getStepSize();
final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
if (nextIsLast) {
setStepSize(finalTime - getStepStart().getTime());
}
}
} while (!isLastStep());
final ODEStateAndDerivative finalState = getStepStart();
setStepStart(null);
setStepSize(Double.NaN);
return finalState;
}
}