Adams-Bashforth methods (in fact due to Adams alone) are explicit * multistep ODE solvers. This implementation is a variation of the classical * one: it uses adaptive stepsize to implement error control, whereas * classical implementations are fixed step size. The value of state vector * at step n+1 is a simple combination of the value at step n and of the * derivatives at steps n, n-1, n-2 ... Depending on the number k of previous * steps one wants to use for computing the next value, different formulas * are available:
*A k-steps Adams-Bashforth method is of order k.
* *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 ** *
The definitions above use the classical representation with several previous first * derivatives. Lets define *
* qn = [ s1(n-1) s1(n-2) ... s1(n-(k-1)) ]T ** (we omit the k index in the notation for clarity). With these definitions, * Adams-Bashforth methods can be written: *
Instead of using the classical representation with first derivatives only (yn, * s1(n) and qn), our 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 ** (here again we omit the k index in the notation for clarity) * * *
Taylor series formulas show that for any index offset i, s1(n-i) can be * computed from s1(n), s2(n) ... sk(n), the formula being exact * for degree k polynomials. *
* s1(n-i) = s1(n) + ∑j j (-i)j-1 sj(n) ** The previous formula can be used with several values for i to compute the transform between * classical representation and Nordsieck vector. The transform between rn * and qn resulting from the Taylor series formulas above is: *
* qn = s1(n) u + P rn ** where u is the [ 1 1 ... 1 ]T vector and P is the (k-1)×(k-1) matrix built * with the j (-i)j-1 terms: *
* [ -2 3 -4 5 ... ] * [ -4 12 -32 80 ... ] * P = [ -6 27 -108 405 ... ] * [ -8 48 -256 1280 ... ] * [ ... ] ** *
Using the Nordsieck vector has several advantages: *
The Nordsieck vector at step n+1 is computed from the Nordsieck vector at step n as follows: *
* [ 0 0 ... 0 0 | 0 ] * [ ---------------+---] * [ 1 0 ... 0 0 | 0 ] * A = [ 0 1 ... 0 0 | 0 ] * [ ... | 0 ] * [ 0 0 ... 1 0 | 0 ] * [ 0 0 ... 0 1 | 0 ] ** *
The P-1u vector and the P-1 A P matrix do not depend on the state, * they only depend on k and therefore are precomputed once for all.
* * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ * @since 2.0 */ public class AdamsBashforthIntegrator extends AdamsIntegrator { /** Integrator method name. */ private static final String METHOD_NAME = "Adams-Bashforth"; /** * Build an Adams-Bashforth integrator with the given order and step control parameters. * @param nSteps number of steps of the method excluding the one being computed * @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 * @exception IllegalArgumentException if order is 1 or less */ public AdamsBashforthIntegrator(final int nSteps, final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) throws IllegalArgumentException { super(METHOD_NAME, nSteps, nSteps, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); } /** * Build an Adams-Bashforth integrator with the given order and step control parameters. * @param nSteps number of steps of the method excluding the one being computed * @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 * @exception IllegalArgumentException if order is 1 or less */ public AdamsBashforthIntegrator(final int nSteps, final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) throws IllegalArgumentException { super(METHOD_NAME, nSteps, nSteps, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); } /** {@inheritDoc} */ @Override public double integrate(final FirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws DerivativeException, IntegratorException { final int n = y0.length; sanityChecks(equations, t0, y0, t, y); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // initialize working arrays if (y != y0) { System.arraycopy(y0, 0, y, 0, n); } final double[] yDot = new double[n]; // set up an interpolator sharing the integrator arrays final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator(); interpolator.reinitialize(y, forward); // set up integration control objects for (StepHandler handler : stepHandlers) { handler.reset(); } setStateInitialized(false); // compute the initial Nordsieck vector using the configured starter integrator start(t0, y, t); interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck); interpolator.storeTime(stepStart); final int lastRow = nordsieck.getRowDimension() - 1; // reuse the step that was chosen by the starter integrator double hNew = stepSize; interpolator.rescale(hNew); // main integration loop isLastStep = false; do { double error = 10; while (error >= 1.0) { stepSize = hNew; // evaluate error using the last term of the Taylor expansion error = 0; for (int i = 0; i < mainSetDimension; ++i) { final double yScale = FastMath.abs(y[i]); final double tol = (vecAbsoluteTolerance == null) ? (scalAbsoluteTolerance + scalRelativeTolerance * yScale) : (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale); final double ratio = nordsieck.getEntry(lastRow, i) / tol; error += ratio * ratio; } error = FastMath.sqrt(error / mainSetDimension); if (error >= 1.0) { // reject the step and attempt to reduce error by stepsize control final double factor = computeStepGrowShrinkFactor(error); hNew = filterStep(stepSize * factor, forward, false); interpolator.rescale(hNew); } } // predict a first estimate of the state at step end final double stepEnd = stepStart + stepSize; interpolator.shift(); interpolator.setInterpolatedTime(stepEnd); System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length); // evaluate the derivative computeDerivatives(stepEnd, y, yDot); // update Nordsieck vector final double[] predictedScaled = new double[y0.length]; for (int j = 0; j < y0.length; ++j) { predictedScaled[j] = stepSize * yDot[j]; } final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck); updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp); interpolator.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp); // discrete events handling interpolator.storeTime(stepEnd); stepStart = acceptStep(interpolator, y, yDot, t); scaled = predictedScaled; nordsieck = nordsieckTmp; 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); final double stopTime = stepStart; resetInternalState(); return stopTime; } }