in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ode/nonstiff/AdamsMoultonFieldIntegrator.java [214:331]
public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
final FieldODEState<T> initialState,
final T finalTime)
throws NumberIsTooSmallException, DimensionMismatchException,
MaxCountExceededException, NoBracketingException {
sanityChecks(initialState, finalTime);
final T t0 = initialState.getTime();
final T[] y = equations.getMapper().mapState(initialState);
setStepStart(initIntegration(equations, t0, y, finalTime));
final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
// compute the initial Nordsieck vector using the configured starter integrator
start(equations, getStepStart(), finalTime);
// reuse the step that was chosen by the starter integrator
FieldODEStateAndDerivative<T> stepStart = getStepStart();
FieldODEStateAndDerivative<T> stepEnd =
AdamsFieldStepInterpolator.taylor(stepStart,
stepStart.getTime().add(getStepSize()),
getStepSize(), scaled, nordsieck);
// main integration loop
setIsLastStep(false);
do {
T[] predictedY = null;
final T[] predictedScaled = MathArrays.buildArray(getField(), y.length);
Array2DRowFieldMatrix<T> predictedNordsieck = null;
T error = getField().getZero().add(10);
while (error.subtract(1.0).getReal() >= 0.0) {
// predict a first estimate of the state at step end (P in the PECE sequence)
predictedY = stepEnd.getState();
// evaluate a first estimate of the derivative (first E in the PECE sequence)
final T[] yDot = computeDerivatives(stepEnd.getTime(), predictedY);
// update Nordsieck vector
for (int j = 0; j < predictedScaled.length; ++j) {
predictedScaled[j] = getStepSize().multiply(yDot[j]);
}
predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);
// apply correction (C in the PECE sequence)
error = predictedNordsieck.walkInOptimizedOrder(new Corrector(y, predictedScaled, predictedY));
if (error.subtract(1.0).getReal() >= 0.0) {
// reject the step and attempt to reduce error by stepsize control
final T factor = computeStepGrowShrinkFactor(error);
rescale(filterStep(getStepSize().multiply(factor), forward, false));
stepEnd = AdamsFieldStepInterpolator.taylor(getStepStart(),
getStepStart().getTime().add(getStepSize()),
getStepSize(),
scaled,
nordsieck);
}
}
// evaluate a final estimate of the derivative (second E in the PECE sequence)
final T[] correctedYDot = computeDerivatives(stepEnd.getTime(), predictedY);
// update Nordsieck vector
final T[] correctedScaled = MathArrays.buildArray(getField(), y.length);
for (int j = 0; j < correctedScaled.length; ++j) {
correctedScaled[j] = getStepSize().multiply(correctedYDot[j]);
}
updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, predictedNordsieck);
// discrete events handling
stepEnd = new FieldODEStateAndDerivative<>(stepEnd.getTime(), predictedY, correctedYDot);
setStepStart(acceptStep(new AdamsFieldStepInterpolator<>(getStepSize(), stepEnd,
correctedScaled, predictedNordsieck, forward,
getStepStart(), stepEnd,
equations.getMapper()),
finalTime));
scaled = correctedScaled;
nordsieck = predictedNordsieck;
if (!isLastStep()) {
System.arraycopy(predictedY, 0, y, 0, y.length);
if (resetOccurred()) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
start(equations, getStepStart(), finalTime);
}
// stepsize control for next step
final T factor = computeStepGrowShrinkFactor(error);
final T scaledH = getStepSize().multiply(factor);
final T nextT = getStepStart().getTime().add(scaledH);
final boolean nextIsLast = forward ?
nextT.subtract(finalTime).getReal() >= 0 :
nextT.subtract(finalTime).getReal() <= 0;
T hNew = filterStep(scaledH, forward, nextIsLast);
final T filteredNextT = getStepStart().getTime().add(hNew);
final boolean filteredNextIsLast = forward ?
filteredNextT.subtract(finalTime).getReal() >= 0 :
filteredNextT.subtract(finalTime).getReal() <= 0;
if (filteredNextIsLast) {
hNew = finalTime.subtract(getStepStart().getTime());
}
rescale(hNew);
stepEnd = AdamsFieldStepInterpolator.taylor(getStepStart(), getStepStart().getTime().add(getStepSize()),
getStepSize(), scaled, nordsieck);
}
} while (!isLastStep());
final FieldODEStateAndDerivative<T> finalState = getStepStart();
setStepStart(null);
setStepSize(null);
return finalState;
}