Added a fast single-step method for fixed-step Runge-Kutta integrators.

JIRA: MATH-1119

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1588769 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2014-04-20 14:29:42 +00:00
parent 2985728c00
commit cbc32459f8
3 changed files with 93 additions and 0 deletions

View File

@ -51,6 +51,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</properties>
<body>
<release version="3.3" date="TBD" description="TBD">
<action dev="luc" type="add" issue="MATH-1119">
Added a fast single-step method for fixed-step Runge-Kutta integrators.
</action>
<action dev="erans" type="fix" issue="MATH-1118">
"Complex": Fixed compatibility of "equals(Object)" with "hashCode()".
Added new methods for testing floating-point equality between the real

View File

@ -24,6 +24,7 @@ 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.FastMath;
/**
@ -185,4 +186,73 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
}
/** 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 = y0.clone();
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 = y0.clone();
// 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;
}
}

View File

@ -306,4 +306,24 @@ public class LutherIntegratorTest {
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
@Test
public void testSingleStep() {
final TestProblem3 pb = new TestProblem3(0.9);
double h = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003;
RungeKuttaIntegrator integ = new LutherIntegrator(Double.NaN);
double t = pb.getInitialTime();
double[] y = pb.getInitialState();
for (int i = 0; i < 100; ++i) {
y = integ.singleStep(pb, t, y, t + h);
t += h;
}
double[] yth = pb.computeTheoreticalState(t);
double dx = y[0] - yth[0];
double dy = y[1] - yth[1];
double error = dx * dx + dy * dy;
Assert.assertEquals(0.0, error, 1.0e-11);
}
}