We define scaled derivatives si(n) at step n as: *
* s1(n) = h y'n for first derivative * s2(n) = h2/2 y''n for second derivative * s3(n) = h3/6 y'''n for third derivative * ... * sk(n) = hk/k! y(k)n for kth derivative **
Rather than storing several previous steps separately, this implementation uses * the Nordsieck vector with higher degrees scaled derivatives all taken at the same * step (yn, s1(n) and rn) where rn is defined as: *
* rn = [ s2(n), s3(n) ... sk(n) ]T ** (we omit the k index in the notation for clarity) *
* Multistep integrators with Nordsieck representation are highly sensitive to * large step changes because when the step is multiplied by a factor a, the * kth component of the Nordsieck vector is multiplied by ak * and the last components are the least accurate ones. The default max growth * factor is therefore set to a quite low value: 21/order. *
* * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ * @since 2.0 */ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { /** First scaled derivative (h y'). */ protected double[] scaled; /** Nordsieck matrix of the higher scaled derivatives. *(h2/2 y'', h3/6 y''' ..., hk/k! y(k))
*/ protected Array2DRowRealMatrix nordsieck; /** Starter integrator. */ private FirstOrderIntegrator starter; /** Number of steps of the multistep method (excluding the one being computed). */ private final int nSteps; /** Stepsize control exponent. */ private double exp; /** Safety factor for stepsize control. */ private double safety; /** Minimal reduction factor for stepsize control. */ private double minReduction; /** Maximal growth factor for stepsize control. */ private double maxGrowth; /** * Build a multistep integrator with the given stepsize bounds. *The default starter integrator is set to the {@link * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with * some defaults settings.
** The default max growth factor is set to a quite low value: 21/order. *
* @param name name of the method * @param nSteps number of steps of the multistep method * (excluding the one being computed) * @param order order of the method * @param minStep minimal step (must be positive even for backward * integration), the last step can be smaller than this * @param maxStep maximal step (must be positive even for backward * integration) * @param scalAbsoluteTolerance allowed absolute error * @param scalRelativeTolerance allowed relative error */ protected MultistepIntegrator(final String name, final int nSteps, final int order, final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); if (nSteps <= 0) { throw MathRuntimeException.createIllegalArgumentException( LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT, name); } starter = new DormandPrince853Integrator(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); this.nSteps = nSteps; exp = -1.0 / order; // set the default values of the algorithm control parameters setSafety(0.9); setMinReduction(0.2); setMaxGrowth(FastMath.pow(2.0, -exp)); } /** * Build a multistep integrator with the given stepsize bounds. *The default starter integrator is set to the {@link * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with * some defaults settings.
** The default max growth factor is set to a quite low value: 21/order. *
* @param name name of the method * @param nSteps number of steps of the multistep method * (excluding the one being computed) * @param order order of the method * @param minStep minimal step (must be positive even for backward * integration), the last step can be smaller than this * @param maxStep maximal step (must be positive even for backward * integration) * @param vecAbsoluteTolerance allowed absolute error * @param vecRelativeTolerance allowed relative error */ protected MultistepIntegrator(final String name, final int nSteps, final int order, final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); starter = new DormandPrince853Integrator(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); this.nSteps = nSteps; exp = -1.0 / order; // set the default values of the algorithm control parameters setSafety(0.9); setMinReduction(0.2); setMaxGrowth(FastMath.pow(2.0, -exp)); } /** * Get the starter integrator. * @return starter integrator */ public ODEIntegrator getStarterIntegrator() { return starter; } /** * Set the starter integrator. *The various step and event handlers for this starter integrator * will be managed automatically by the multi-step integrator. Any * user configuration for these elements will be cleared before use.
* @param starterIntegrator starter integrator */ public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) { this.starter = starterIntegrator; } /** Start the integration. *This method computes one step using the underlying starter integrator, * and initializes the Nordsieck vector at step start. The starter integrator * purpose is only to establish initial conditions, it does not really change * time by itself. The top level multistep integrator remains in charge of * handling time propagation and events handling as it will starts its own * computation right from the beginning. In a sense, the starter integrator * can be seen as a dummy one and so it will never trigger any user event nor * call any user step handler.
* @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 thant0 for backward integration)
* @throws IntegratorException if the integrator cannot perform integration
* @throws DerivativeException this exception is propagated to the caller if
* the underlying user function triggers one
*/
protected void start(final double t0, final double[] y0, final double t)
throws DerivativeException, IntegratorException {
// make sure NO user event nor user step handler is triggered,
// this is the task of the top level integrator, not the task
// of the starter integrator
starter.clearEventHandlers();
starter.clearStepHandlers();
// set up one specific step handler to extract initial Nordsieck vector
starter.addStepHandler(new NordsieckInitializer(y0.length));
// start integration, expecting a InitializationCompletedMarkerException
try {
starter.integrate(new CountingDifferentialEquations(y0.length),
t0, y0, t, new double[y0.length]);
} catch (DerivativeException mue) {
if (!(mue instanceof InitializationCompletedMarkerException)) {
// this is not the expected nominal interruption of the start integrator
throw mue;
}
}
// remove the specific step handler
starter.clearStepHandlers();
}
/** Initialize the high order scaled derivatives at step start.
* @param first first scaled derivative at step start
* @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
* will be modified
* @return high order scaled derivatives at step start
*/
protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
final double[][] multistep);
/** Get the minimal reduction factor for stepsize control.
* @return minimal reduction factor
*/
public double getMinReduction() {
return minReduction;
}
/** Set the minimal reduction factor for stepsize control.
* @param minReduction minimal reduction factor
*/
public void setMinReduction(final double minReduction) {
this.minReduction = minReduction;
}
/** Get the maximal growth factor for stepsize control.
* @return maximal growth factor
*/
public double getMaxGrowth() {
return maxGrowth;
}
/** Set the maximal growth factor for stepsize control.
* @param maxGrowth maximal growth factor
*/
public void setMaxGrowth(final double maxGrowth) {
this.maxGrowth = maxGrowth;
}
/** Get the safety factor for stepsize control.
* @return safety factor
*/
public double getSafety() {
return safety;
}
/** Set the safety factor for stepsize control.
* @param safety safety factor
*/
public void setSafety(final double safety) {
this.safety = safety;
}
/** Compute step grow/shrink factor according to normalized error.
* @param error normalized error of the current step
* @return grow/shrink factor for next step
*/
protected double computeStepGrowShrinkFactor(final double error) {
return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
}
/** Transformer used to convert the first step to Nordsieck representation. */
public static interface NordsieckTransformer {
/** Initialize the high order scaled derivatives at step start.
* @param first first scaled derivative at step start
* @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
* will be modified
* @return high order derivatives at step start
*/
RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
}
/** Specialized step handler storing the first step. */
private class NordsieckInitializer implements StepHandler {
/** Problem dimension. */
private final int n;
/** Simple constructor.
* @param n problem dimension
*/
public NordsieckInitializer(final int n) {
this.n = n;
}
/** {@inheritDoc} */
public void handleStep(StepInterpolator interpolator, boolean isLast)
throws DerivativeException {
final double prev = interpolator.getPreviousTime();
final double curr = interpolator.getCurrentTime();
stepStart = prev;
stepSize = (curr - prev) / (nSteps + 1);
// compute the first scaled derivative
interpolator.setInterpolatedTime(prev);
scaled = interpolator.getInterpolatedDerivatives().clone();
for (int j = 0; j < n; ++j) {
scaled[j] *= stepSize;
}
// compute the high order scaled derivatives
final double[][] multistep = new double[nSteps][];
for (int i = 1; i <= nSteps; ++i) {
interpolator.setInterpolatedTime(prev + stepSize * i);
final double[] msI = interpolator.getInterpolatedDerivatives().clone();
for (int j = 0; j < n; ++j) {
msI[j] *= stepSize;
}
multistep[i - 1] = msI;
}
nordsieck = initializeHighOrderDerivatives(scaled, multistep);
// stop the integrator after the first step has been handled
throw new InitializationCompletedMarkerException();
}
/** {@inheritDoc} */
public boolean requiresDenseOutput() {
return true;
}
/** {@inheritDoc} */
public void reset() {
// nothing to do
}
}
/** Marker exception used ONLY to stop the starter integrator after first step. */
private static class InitializationCompletedMarkerException
extends DerivativeException {
/** Serializable version identifier. */
private static final long serialVersionUID = -4105805787353488365L;
/** Simple constructor. */
public InitializationCompletedMarkerException() {
super((Throwable) null);
}
}
/** Wrapper for differential equations, ensuring start evaluations are counted. */
private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations {
/** Dimension of the problem. */
private final int dimension;
/** Simple constructor.
* @param dimension dimension of the problem
*/
public CountingDifferentialEquations(final int dimension) {
this.dimension = dimension;
}
/** {@inheritDoc} */
public void computeDerivatives(double t, double[] y, double[] dot)
throws DerivativeException {
MultistepIntegrator.this.computeDerivatives(t, y, dot);
}
/** {@inheritDoc} */
public int getDimension() {
return dimension;
}
/** {@inheritDoc} */
public int getMainSetDimension() {
return mainSetDimension;
}
}
}