in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ode/nonstiff/AdamsBashforthIntegrator.java [234:353]
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();
// initialize working arrays
final double[] y = equations.getCompleteState();
final double[] yDot = new double[y.length];
// set up an interpolator sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
// set up integration control objects
initIntegration(equations.getTime(), y, t);
// compute the initial Nordsieck vector using the configured starter integrator
start(equations.getTime(), y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
// reuse the step that was chosen by the starter integrator
double hNew = stepSize;
interpolator.rescale(hNew);
// main integration loop
isLastStep = false;
do {
interpolator.shift();
final double[] predictedY = new double[y.length];
final double[] predictedScaled = new double[y.length];
Array2DRowRealMatrix predictedNordsieck = null;
double error = 10;
while (error >= 1.0) {
// predict a first estimate of the state at step end
final double stepEnd = stepStart + hNew;
interpolator.storeTime(stepEnd);
final ExpandableStatefulODE expandable = getExpandable();
final EquationsMapper primary = expandable.getPrimaryMapper();
primary.insertEquationData(interpolator.getInterpolatedState(), predictedY);
int index = 0;
for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), predictedY);
++index;
}
// evaluate the derivative
computeDerivatives(stepEnd, predictedY, yDot);
// predict Nordsieck vector at step end
for (int j = 0; j < predictedScaled.length; ++j) {
predictedScaled[j] = hNew * yDot[j];
}
predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);
// evaluate error
error = errorEstimation(y, predictedY, predictedScaled, predictedNordsieck);
if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor = computeStepGrowShrinkFactor(error);
hNew = filterStep(hNew * factor, forward, false);
interpolator.rescale(hNew);
}
}
stepSize = hNew;
final double stepEnd = stepStart + stepSize;
interpolator.reinitialize(stepEnd, stepSize, predictedScaled, predictedNordsieck);
// discrete events handling
interpolator.storeTime(stepEnd);
System.arraycopy(predictedY, 0, y, 0, y.length);
stepStart = acceptStep(interpolator, y, yDot, t);
scaled = predictedScaled;
nordsieck = predictedNordsieck;
interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
if (!isLastStep) {
// prepare next step
interpolator.storeTime(stepStart);
if (resetOccurred) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
start(stepStart, y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
}
// stepsize control for next step
final double factor = computeStepGrowShrinkFactor(error);
final double scaledH = stepSize * factor;
final double nextT = stepStart + scaledH;
final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
hNew = filterStep(scaledH, forward, nextIsLast);
final double filteredNextT = stepStart + hNew;
final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
if (filteredNextIsLast) {
hNew = t - stepStart;
}
interpolator.rescale(hNew);
}
} while (!isLastStep);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
resetInternalState();
}