commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r676737 [1/2] - in /commons/proper/math/branches/MATH_2_0/src: java/org/apache/commons/math/ode/nonstiff/ site/xdoc/ test/org/apache/commons/math/ode/nonstiff/
Date Mon, 14 Jul 2008 21:13:23 GMT
Author: luc
Date: Mon Jul 14 14:13:23 2008
New Revision: 676737

URL: http://svn.apache.org/viewvc?rev=676737&view=rev
Log:
New ODE integrators have been added: the explicit Adams-Bashforth and implicit
Adams-Moulton multistep methods. These methods support customizable starter
integrators and support discrete events even during the start phase.
All these methods provide the same rich features has the existing ones:
continuous output, step handlers, discrete events, G-stop ...

Added:
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java   (with props)
    commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java   (with props)
Modified:
    commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,285 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.sampling.StepHandler;
+
+
+/**
+ * This class implements explicit Adams-Bashforth integrators for Ordinary
+ * Differential Equations.
+ *
+ * <p>Adams-Bashforth (in fact due to Adams alone) methods are explicit
+ * multistep ODE solvers witch fixed stepsize. 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:</p>
+ * <ul>
+ *   <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n</sub></li>
+ *   <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (3f<sub>n</sub>-f<sub>n-1</sub>)/2</li>
+ *   <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (23f<sub>n</sub>-16f<sub>n-1</sub>+5f<sub>n-2</sub>)/12</li>
+ *   <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + h (55f<sub>n</sub>-59f<sub>n-1</sub>+37f<sub>n-2</sub>-9f<sub>n-3)/24</sub></li>
+ *   <li>...</li>
+ * </ul>
+ *
+ * <p>A k-steps Adams-Bashforth method is of order k. There is no limit to the
+ * value of k.</p>
+ *
+ * <p>These methods are explicit: f<sub>n+1</sub> is not used to compute
+ * y<sub>n+1</sub>. More accurate implicit Adams methods exist: the
+ * Adams-Moulton methods (which are also due to Adams alone). They are
+ * provided by the {@link AdamsMoultonIntegrator AdamsMoultonIntegrator} class.</p>
+ *
+ * @see AdamsMoultonIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class AdamsBashforthIntegrator extends MultistepIntegrator {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 1676381657635800870L;
+
+    /** Integrator method name. */
+    private static final String METHOD_NAME = "Adams-Bashforth";
+
+   /** Coefficients for the current method. */
+    private final double[] coeffs;
+
+    /** Integration step. */
+    private final double step;
+
+    /**
+     * Build an Adams-Bashforth integrator with the given order and step size.
+     * @param order order of the method (must be strictly positive)
+     * @param step integration step size
+     */
+    public AdamsBashforthIntegrator(final int order, final double step) {
+
+        super(METHOD_NAME, order, new AdamsBashforthStepInterpolator());
+
+        // compute the integration coefficients
+        int[][] bdArray = computeBackwardDifferencesArray(order);
+        Fraction[] gamma = computeGammaArray(order);
+        coeffs = new double[order];
+        for (int i = 0; i < order; ++i) {
+            Fraction f = Fraction.ZERO;
+            for (int j = i; j < order; ++j) {
+                f = f.add(gamma[j].multiply(new Fraction(bdArray[j][i], 1)));
+            }
+            coeffs[i] = f.doubleValue();
+        }
+
+        this.step = step;
+
+    }
+
+    /** {@inheritDoc} */
+    public double integrate(FirstOrderDifferentialEquations equations,
+                            double t0, double[] y0, double t, double[] y)
+        throws DerivativeException, IntegratorException {
+
+        sanityChecks(equations, t0, y0, t, y);
+        final boolean forward = (t > t0);
+
+        // initialize working arrays
+        if (y != y0) {
+            System.arraycopy(y0, 0, y, 0, y0.length);
+        }
+        final double[] yTmp = new double[y0.length];
+
+        // set up an interpolator sharing the integrator arrays
+        final AdamsBashforthStepInterpolator interpolator =
+                (AdamsBashforthStepInterpolator) prototype.copy();
+        interpolator.reinitialize(yTmp, previousT, previousF, forward);
+
+        // set up integration control objects
+        stepStart = t0;
+        stepSize  = step;
+        for (StepHandler handler : stepHandlers) {
+            handler.reset();
+        }
+        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+
+        // compute the first few steps using the configured starter integrator
+        double stopTime =
+            start(previousF.length, stepSize, manager, equations, stepStart, y);
+        if (Double.isNaN(previousT[0])) {
+            return stopTime;
+        }
+        stepStart = previousT[0];
+        interpolator.storeTime(stepStart);
+
+        boolean lastStep = false;
+        while (!lastStep) {
+
+            // shift all data
+            interpolator.shift();
+
+            // estimate the state at the end of the step
+            for (int j = 0; j < y0.length; ++j) {
+                double sum = 0;
+                for (int l = 0; l < coeffs.length; ++l) {
+                    sum += coeffs[l] * previousF[l][j];
+                }
+                yTmp[j] = y[j] + stepSize * sum;
+            }
+
+            // discrete events handling
+            interpolator.storeTime(stepStart + stepSize);
+            final boolean truncated;
+            if (manager.evaluateStep(interpolator)) {
+                truncated = true;
+                interpolator.truncateStep(manager.getEventTime());
+            } else {
+                truncated = false;
+            }
+
+            // the step has been accepted (may have been truncated)
+            final double nextStep = interpolator.getCurrentTime();
+            interpolator.setInterpolatedTime(nextStep);
+            System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+            manager.stepAccepted(nextStep, y);
+            lastStep = manager.stop();
+
+            // provide the step data to the step handler
+            for (StepHandler handler : stepHandlers) {
+                handler.handleStep(interpolator, lastStep);
+            }
+            stepStart = nextStep;
+
+            if (!lastStep) {
+                // prepare next step
+
+                if (manager.reset(stepStart, y)) {
+
+                    // some events handler has triggered changes that
+                    // invalidate the derivatives, we need to restart from scratch
+                    stopTime =
+                        start(previousF.length, stepSize, manager, equations, stepStart, y);
+                    if (Double.isNaN(previousT[0])) {
+                        return stopTime;
+                    }
+                    stepStart = previousT[0];
+
+                } else {
+
+                    if (truncated) {
+                        // the step has been truncated, we need to adjust the previous steps
+                        for (int i = 1; i < previousF.length; ++i) {
+                            previousT[i] = stepStart - i * stepSize;
+                            interpolator.setInterpolatedTime(previousT[i]);
+                            System.arraycopy(interpolator.getInterpolatedState(), 0,
+                                             previousF[i], 0, y0.length);
+                        }
+                    } else {
+                        rotatePreviousSteps();
+                    }
+
+                    // evaluate differential equations for next step
+                    previousT[0] = stepStart;
+                    equations.computeDerivatives(stepStart, y, previousF[0]);
+
+                }
+            }
+
+        }
+
+        stopTime  = stepStart;
+        stepStart = Double.NaN;
+        stepSize  = Double.NaN;
+        return stopTime;
+
+    }
+
+    /** Get the coefficients of the method.
+     * <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
+     * <pre>
+     *   y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k-1</sup> c<sub>i</sub>f<sub>n-i</sub></li>
+     * </pre>
+     * @return a copy of the coefficients of the method
+     */
+    public double[] getCoeffs() {
+        return coeffs.clone();
+    }
+
+    /** Compute the backward differences coefficients array.
+     * <p>This is quite similar to the Pascal triangle, except for a
+     * (-1)<sup>i</sup> sign. We use a straightforward approach here,
+     * since we don't expect this to be run too many times with too
+     * high k. It is based on the recurrence relations:</p>
+     * <pre>
+     *   &nabla;<sup>0</sup> f<sub>n</sub> = f<sub>n</sub>
+     *   &nabla;<sup>i+1</sup> f<sub>n</sub> = &nabla;<sup>i</sup>f<sub>n</sub> - &nabla;<sup>i</sup>f<sub>n-1</sub>
+     * </pre>
+     * @param order order of the integration method
+     */
+    static int[][] computeBackwardDifferencesArray(final int order) {
+
+        // create the array
+        int[][] bdArray = new int[order][];
+
+        // recurrence initialization
+        bdArray[0] = new int[] { 1 };
+
+        // fill up array using recurrence relation
+        for (int i = 1; i < order; ++i) {
+            bdArray[i] = new int[i + 1];
+            bdArray[i][0] = 1;
+            for (int j = 0; j < i - 1; ++j) {
+                bdArray[i][j + 1] = bdArray[i - 1][j + 1] - bdArray[i - 1][j];
+            }
+            bdArray[i][i] = -bdArray[i - 1][i - 1];
+        }
+
+        return bdArray;
+
+    }
+
+    /** Compute the gamma coefficients.
+     * @param order order of the integration method
+     * @return gamma coefficients array
+     */
+    static Fraction[] computeGammaArray(final int order) {
+
+        // create the array
+        Fraction[] gammaArray = new Fraction[order];
+
+        // recurrence initialization
+        gammaArray[0] = Fraction.ONE;
+
+        // fill up array using recurrence relation
+        for (int i = 1; i < order; ++i) {
+            Fraction gamma = Fraction.ONE;
+            for (int j = 1; j <= i; ++j) {
+                gamma = gamma.subtract(gammaArray[i - j].multiply(new Fraction(1, j + 1)));
+            }
+            gammaArray[i] = gamma;
+        }
+
+        return gammaArray;
+
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,288 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+
+/**
+ * This class implements an interpolator for Adams-Bashforth multiple steps.
+ *
+ * <p>This interpolator computes dense output inside the last few
+ * steps computed. The interpolation equation is consistent with the
+ * integration scheme, it is based on a kind of <em>rollback</em> of the
+ * integration from step end to interpolation date:
+ * <pre>
+ *   y(t<sub>n</sub> + theta h) = y (t<sub>n</sub> + h) - &int;<sub>t<sub>n</sub> + theta h</sub><sup>t<sub>n</sub> + h</sup>p(t)dt
+ * </pre>
+ * where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on
+ * the derivatives at previous steps f<sub>n-k+1</sub>, f<sub>n-k+2</sub> ...
+ * f<sub>n</sub> and f<sub>n</sub>.</p>
+ *
+ * @see AdamsBashforthIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+class AdamsBashforthStepInterpolator extends MultistepStepInterpolator {
+
+    /** Serializable version identifier */
+    private static final long serialVersionUID = -7179861704951334960L;
+
+    /** Neville's interpolation array. */
+    private double[] neville;
+
+    /** Integration rollback array. */
+    private double[] rollback;
+
+    /** &gamma; array. */
+    private double[] gamma;
+
+    /** Backward differences array. */
+    private int[][] bdArray;
+
+    /** Original non-truncated step end time. */
+    private double nonTruncatedEnd;
+
+    /** Original non-truncated step size. */
+    private double nonTruncatedH;
+
+    /** Simple constructor.
+     * This constructor builds an instance that is not usable yet, the
+     * {@link AbstractStepInterpolator#reinitialize} method should be called
+     * before using the instance in order to initialize the internal arrays. This
+     * constructor is used only in order to delay the initialization in
+     * some cases.
+     */
+    public AdamsBashforthStepInterpolator() {
+    }
+
+    /** Copy constructor.
+     * @param interpolator interpolator to copy from. The copy is a deep
+     * copy: its arrays are separated from the original arrays of the
+     * instance
+     */
+    public AdamsBashforthStepInterpolator(final AdamsBashforthStepInterpolator interpolator) {
+        super(interpolator);
+        nonTruncatedEnd = interpolator.nonTruncatedEnd;
+        nonTruncatedH   = interpolator.nonTruncatedH;
+    }
+
+    /** {@inheritDoc} */
+    protected StepInterpolator doCopy() {
+        return new AdamsBashforthStepInterpolator(this);
+    }
+
+    /** {@inheritDoc} */
+    protected void initializeCoefficients() {
+
+        neville  = new double[previousF.length];
+        rollback = new double[previousF.length];
+
+        bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
+
+        Fraction[] fGamma = AdamsBashforthIntegrator.computeGammaArray(previousF.length);
+        gamma = new double[fGamma.length];
+        for (int i = 0; i < fGamma.length; ++i) {
+            gamma[i] = fGamma[i].doubleValue();
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public void storeTime(final double t) {
+        nonTruncatedEnd = t;
+        nonTruncatedH   = nonTruncatedEnd - previousTime;
+        super.storeTime(t);
+    }
+
+    /** Truncate a step.
+     * <p>Truncating a step is necessary when an event is triggered
+     * before the nominal end of the step.</p>
+     */
+    void truncateStep(final double truncatedEndTime) {
+        currentTime = truncatedEndTime;
+        h = currentTime - previousTime;
+    }
+
+    /** {@inheritDoc} */
+    public void setInterpolatedTime(final double time)
+        throws DerivativeException {
+        interpolatedTime = time;
+        final double oneMinusThetaH = nonTruncatedEnd - interpolatedTime;
+        final double theta = (nonTruncatedH == 0) ?
+                             0 : (nonTruncatedH - oneMinusThetaH) / nonTruncatedH;
+        computeInterpolatedState(theta, oneMinusThetaH);
+    }
+
+    /** {@inheritDoc} */
+    protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
+        interpolateDerivatives();
+        interpolateState(theta);
+    }
+
+    /** Interpolate the derivatives.
+     * <p>The Adams method is based on a polynomial interpolation of the
+     * derivatives based on the preceding steps. So the interpolation of
+     * the derivatives here is strictly equivalent: it is a simple polynomial
+     * interpolation.</p>
+     */
+    private void interpolateDerivatives() {
+
+        for (int i = 0; i < interpolatedDerivatives.length; ++i) {
+
+            // initialize the Neville's interpolation algorithm
+            for (int k = 0; k < previousF.length; ++k) {
+                neville[k] = previousF[k][i];
+            }
+
+            // combine the contributions of each points
+            for (int l = 1; l < neville.length; ++l) {
+                for (int m = neville.length - 1; m >= l; --m) {
+                    final double xm   = previousT[m];
+                    final double xmMl = previousT[m - l];
+                    neville[m] = ((interpolatedTime - xm) * neville[m-1] +
+                                  (xmMl - interpolatedTime) * neville[m]) / (xmMl - xm);
+                }
+            }
+
+            // the interpolation polynomial value is in the array last element
+            interpolatedDerivatives[i] = neville[neville.length - 1];
+
+        }
+
+    }
+
+    /** Interpolate the state.
+     * <p>The Adams method is based on a polynomial interpolation of the
+     * derivatives based on the preceding steps. The polynomial model is
+     * integrated analytically throughout the last step. Using the notations
+     * found in the second edition of the first volume (Nonstiff Problems)
+     * of the reference book by Hairer, Norsett and Wanner: <i>Solving Ordinary
+     * Differential Equations</i> (Springer-Verlag, ISBN 3-540-56670-8), this
+     * process leads to the following expression:</p>
+     * <pre>
+     * y<sub>n+1</sub> = y<sub>n</sub> +
+     * h &times; &sum;<sub>j=0</sub><sup>j=k-1</sup> &gamma;<sub>j</sub>&nabla;<sup>j</sup>f<sub>n</sub>
+     * </pre>
+     * <p>In the previous expression, the &gamma;<sub>j</sub> terms are the
+     * ones that result from the analytical integration, and can be computed form
+     * the binomial coefficients C<sub>j</sub><sup>-s</sup>:</p>
+     * <p>
+     * &gamma;<sub>j</sub> = (-1)<sup>j</sup>&int;<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
+     * </p>
+     * <p>In order to interpolate the state in a manner that is consistent with the
+     * integration scheme, we simply subtract from the current state (at the end of the step)
+     * the integral computed from interpolation time to step end time.</p>
+     * <p>
+     * &eta;<sub>j</sub>(&theta;)=
+     * (-1)<sup>j</sup>&int;<sub>&theta;</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
+     * </p>
+     * The method described in the Hairer, Norsett and Wanner book to compute &gamma;<sub>j</sub>
+     * is easily extended to compute &gamma;<sub>j</sub>(&theta;)=
+     * (-1)<sup>j</sup>&int;<sub>0</sub><sup>&theta;</sup>C<sub>j</sub><sup>-s</sup>ds. From this,
+     * we can compute &eta;<sub>j</sub>(&theta;) = &gamma;<sub>j</sub>-&gamma;<sub>j</sub>(&theta;).
+     * The first few values are:</p>
+     * <table>
+     * <tr><td>j</td><td>&gamma;<sub>j</sub></td><td>&gamma;<sub>j</sub>(&theta;)</td><td>&eta;<sub>j</sub>(&theta;)</td></tr>
+     * <tr><td>0</td><td>1</td><td></td>&theta;<td>1-&theta;</td></tr>
+     * <tr><td>1</td><td>1/2</td><td></td>&theta;<sup>2</sup>/2<td>(1-&theta;<sup>2</sup>)/2</td></tr>
+     * <tr><td>2</td><td>5/12</td><td></td>(3&theta;<sup>2</sup>+2&theta;<sup>3</sup>)/12<td>(5-3&theta;<sup>2</sup>-2&theta;<sup>3</sup>)/12</td></tr>
+     * </table>
+     * <p>
+     * The &eta;<sub>j</sub>(&theta;) functions appear to be polynomial ones. As expected,
+     * we see that &eta;<sub>j</sub>(1)= 0. The recurrence relation derived for
+     * &gamma;<sub>j</sub>(&theta;) is:
+     * </p>
+     * <p>
+     * &sum<sub>j=0</sub><sup>j=m</sup>&gamma;<sub>j</sub>(&theta;)/(m+1-j) =
+     * 1/(m+1)! &prod;<sub>k=0</sub><sup>k=m</sup>(&theta;+k)
+     * </p>
+     * @param theta location of the interpolation point within the last step
+     */
+    private void interpolateState(final double theta) {
+
+        // compute the integrals to remove from the final state
+        computeRollback(previousT.length - 1, theta);
+
+        // remove these integrals from the final state
+        for (int j = 0; j < interpolatedState.length; ++j) {
+            double sum = 0;
+            for (int l = 0; l < previousT.length; ++l) {
+                sum += rollback[l] * previousF[l][j];
+            }
+            interpolatedState[j] = currentState[j] - h * sum;
+        }
+
+    }
+
+    /** Compute the rollback coefficients.
+     * @param order order of the integration method
+     * @param theta current value for theta
+     */
+    private void computeRollback(final int order, final double theta) {
+
+        // compute the gamma(theta) values from the recurrence relation
+        double product = theta;
+        rollback[0]  = theta;
+        for (int i = 1; i < order; ++i) {
+            product *= (i + theta) / (i + 1);
+            double g = product;
+            for (int j = 1; j <= i; ++j) {
+                g -= rollback[i - j] / (j + 1);
+            }
+            rollback[i] = g;
+        }
+
+        // subtract it from gamma to get eta(theta)
+        for (int i = 0; i < order; ++i) {
+            rollback[i] -= gamma[i];
+        }
+
+        // combine the eta integrals with the backward differences array
+        // to get the rollback coefficients
+        for (int i = 0; i < order; ++i) {
+            double f = 0;
+            for (int j = i; j <= order; ++j) {
+                f -= rollback[j] * bdArray[j][i];
+            }
+            rollback[i] = f;
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public void writeExternal(final ObjectOutput out)
+        throws IOException {
+        super.writeExternal(out);
+        out.writeDouble(nonTruncatedEnd);
+    }
+
+    /** {@inheritDoc} */
+    public void readExternal(final ObjectInput in)
+        throws IOException {
+        nonTruncatedEnd = in.readDouble();
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,299 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.sampling.StepHandler;
+
+
+/**
+ * This class implements implicit Adams-Moulton integrators for Ordinary
+ * Differential Equations.
+ *
+ * <p>Adams-Moulton (in fact due to Adams alone) methods are implicit
+ * multistep ODE solvers witch fixed stepsize. 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+1, n, n-1 ... Depending on the number k of previous
+ * steps one wants to use for computing the next value, different formulas
+ * are available:</p>
+ * <ul>
+ *   <li>k = 0: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n+1</sub></li>
+ *   <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h (f<sub>n+1</sub>+f<sub>n</sub>)/2</li>
+ *   <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (5f<sub>n+1</sub>+8f<sub>n</sub>-f<sub>n-1</sub>)/12</li>
+ *   <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (9f<sub>n+1</sub>+19f<sub>n</sub>-5f<sub>n-1</sub>+f<sub>n-2)/24</sub></li>
+ *   <li>...</li>
+ * </ul>
+ *
+ * <p>The coefficients are computed (and cached) as needed, so their are no
+ * theoretical limitations to the number of steps</p>
+ *
+ * <p>A k-steps Adams-Moulton method is of order k+1. There is no limit to the
+ * value of k.</p>
+ *
+ * <p>These methods are implicit: f<sub>n+1</sub> is used to compute
+ * y<sub>n+1</sub>. Simpler explicit Adams methods exist: the
+ * Adams-Bashforth methods (which are also due to Adams alone). They are
+ * provided by the {@link AdamsBashforthIntegrator AdamsBashforthIntegrator} class.</p>
+ *
+ * @see AdamsBashforthIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class AdamsMoultonIntegrator extends MultistepIntegrator {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 4990335331377040417L;
+
+    /** Integrator method name. */
+    private static final String METHOD_NAME = "Adams-Moulton";
+
+    /** Coefficients for the predictor phase of the method. */
+    private final double[] predictorCoeffs;
+
+    /** Coefficients for the corrector phase of the method. */
+    private final double[] correctorCoeffs;
+
+    /** Integration step. */
+    private final double step;
+
+    /**
+     * Build an Adams-Moulton integrator with the given order and step size.
+     * @param order order of the method (must be strictly positive)
+     * @param step integration step size
+     */
+    public AdamsMoultonIntegrator(final int order, final double step) {
+
+        super(METHOD_NAME, order + 1, new AdamsMoultonStepInterpolator());
+
+        // compute the integration coefficients
+        int[][] bdArray      = AdamsBashforthIntegrator.computeBackwardDifferencesArray(order + 1);
+
+        Fraction[] gamma     = AdamsBashforthIntegrator.computeGammaArray(order);
+        predictorCoeffs = new double[order];
+        for (int i = 0; i < order; ++i) {
+            Fraction fPredictor = Fraction.ZERO;
+            for (int j = i; j < order; ++j) {
+                Fraction f = new Fraction(bdArray[j][i], 1);
+                fPredictor = fPredictor.add(gamma[j].multiply(f));
+            }
+            predictorCoeffs[i] = fPredictor.doubleValue();
+        }
+
+        Fraction[] gammaStar = computeGammaStarArray(order);
+        correctorCoeffs = new double[order + 1];
+        for (int i = 0; i <= order; ++i) {
+            Fraction fCorrector = Fraction.ZERO;
+            for (int j = i; j <= order; ++j) {
+                Fraction f = new Fraction(bdArray[j][i], 1);
+                fCorrector = fCorrector.add(gammaStar[j].multiply(f));
+            }
+            correctorCoeffs[i] = fCorrector.doubleValue();
+        }
+
+        this.step = step;
+
+    }
+
+    /** {@inheritDoc} */
+    public double integrate(FirstOrderDifferentialEquations equations,
+                            double t0, double[] y0, double t, double[] y)
+        throws DerivativeException, IntegratorException {
+
+        sanityChecks(equations, t0, y0, t, y);
+        final boolean forward = (t > t0);
+
+        // initialize working arrays
+        if (y != y0) {
+            System.arraycopy(y0, 0, y, 0, y0.length);
+        }
+        final double[] yTmp = new double[y0.length];
+
+        // set up an interpolator sharing the integrator arrays
+        final AdamsMoultonStepInterpolator interpolator =
+                (AdamsMoultonStepInterpolator) prototype.copy();
+        interpolator.reinitialize(yTmp, previousT, previousF, forward);
+
+        // set up integration control objects
+        stepStart = t0;
+        stepSize  = step;
+        for (StepHandler handler : stepHandlers) {
+            handler.reset();
+        }
+        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+
+        // compute the first few steps using the configured starter integrator
+        double stopTime =
+            start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
+        if (Double.isNaN(previousT[0])) {
+            return stopTime;
+        }
+        stepStart = previousT[0];
+        rotatePreviousSteps();
+        previousF[0] = new double[y0.length];
+        interpolator.storeTime(stepStart);
+
+        boolean lastStep = false;
+        while (!lastStep) {
+
+            // shift all data
+            interpolator.shift();
+
+            // predict state at end of step
+            for (int j = 0; j < y0.length; ++j) {
+                double sum = 0;
+                for (int l = 0; l < predictorCoeffs.length; ++l) {
+                    sum += predictorCoeffs[l] * previousF[l+1][j];
+                }
+                yTmp[j] = y[j] + stepSize * sum;
+            }
+
+            // evaluate the derivatives
+            final double stepEnd = stepStart + stepSize;
+            equations.computeDerivatives(stepEnd, yTmp, previousF[0]);
+
+            // apply corrector
+            for (int j = 0; j < y0.length; ++j) {
+                double sum = 0;
+                for (int l = 0; l < correctorCoeffs.length; ++l) {
+                    sum += correctorCoeffs[l] * previousF[l][j];
+                }
+                yTmp[j] = y[j] + stepSize * sum;
+            }
+
+            // discrete events handling
+            interpolator.storeTime(stepEnd);
+            final boolean truncated;
+            if (manager.evaluateStep(interpolator)) {
+                truncated = true;
+                interpolator.truncateStep(manager.getEventTime());
+            } else {
+                truncated = false;
+            }
+
+            // the step has been accepted (may have been truncated)
+            final double nextStep = interpolator.getCurrentTime();
+            interpolator.setInterpolatedTime(nextStep);
+            System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+            manager.stepAccepted(nextStep, y);
+            lastStep = manager.stop();
+
+            // provide the step data to the step handler
+            for (StepHandler handler : stepHandlers) {
+                handler.handleStep(interpolator, lastStep);
+            }
+            stepStart = nextStep;
+
+            if (!lastStep) {
+                // prepare next step
+
+                if (manager.reset(stepStart, y)) {
+
+                    // some events handler has triggered changes that
+                    // invalidate the derivatives, we need to restart from scratch
+                    stopTime =
+                        start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
+                    if (Double.isNaN(previousT[0])) {
+                        return stopTime;
+                    }
+                    stepStart = previousT[0];
+                    rotatePreviousSteps();
+                    previousF[0] = new double[y0.length];
+
+                } else {
+
+                    if (truncated) {
+                        // the step has been truncated, we need to adjust the previous steps
+                        for (int i = 1; i < previousF.length; ++i) {
+                            previousT[i] = stepStart - i * stepSize;
+                            interpolator.setInterpolatedTime(previousT[i]);
+                            System.arraycopy(interpolator.getInterpolatedState(), 0,
+                                             previousF[i], 0, y0.length);
+                        }
+                    } else {
+                        rotatePreviousSteps();
+                    }
+
+                    // evaluate differential equations for next step
+                    previousT[0] = stepStart;
+                    equations.computeDerivatives(stepStart, y, previousF[0]);
+
+                }
+            }
+
+        }
+
+        stopTime  = stepStart;
+        stepStart = Double.NaN;
+        stepSize  = Double.NaN;
+        return stopTime;
+
+    }
+
+    /** Get the coefficients of the prdictor phase of the method.
+     * <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
+     * <pre>
+     *   y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k-1</sup> c<sub>i</sub>f<sub>n-i</sub></li>
+     * </pre>
+     * @return a copy of the coefficients of the method
+     */
+    public double[] getPredictorCoeffs() {
+        return predictorCoeffs.clone();
+    }
+
+    /** Get the coefficients of the corrector phase of the method.
+     * <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
+     * <pre>
+     *   y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k</sup> c<sub>i</sub>f<sub>n-i</sub></li>
+     * </pre>
+     * @return a copy of the coefficients of the method
+     */
+    public double[] getCorrectorCoeffs() {
+        return correctorCoeffs.clone();
+    }
+
+    /** Compute the gamma star coefficients.
+     * @param order order of the integration method
+     * @return gamma star coefficients array
+     */
+    static Fraction[] computeGammaStarArray(final int order) {
+
+        // create the array
+        Fraction[] gammaStarArray = new Fraction[order + 1];
+
+        // recurrence initialization
+        gammaStarArray[0] = Fraction.ONE;
+
+        // fill up array using recurrence relation
+        for (int i = 1; i <= order; ++i) {
+            Fraction gammaStar = Fraction.ZERO;
+            for (int j = 1; j <= i; ++j) {
+                gammaStar = gammaStar.subtract(gammaStarArray[i - j].multiply(new Fraction(1, j + 1)));
+            }
+            gammaStarArray[i] = gammaStar;
+        }
+
+        return gammaStarArray;
+
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,289 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+
+/**
+ * This class implements an interpolator for Adams-Moulton multiple steps.
+ *
+ * <p>This interpolator computes dense output inside the last few
+ * steps computed. The interpolation equation is consistent with the
+ * integration scheme, it is based on a kind of <em>rollback</em> of the
+ * integration from step end to interpolation date:
+ * <pre>
+ *   y(t<sub>n</sub> + theta h) = y (t<sub>n</sub> + h) - &int;<sub>t<sub>n</sub> + theta h</sub><sup>t<sub>n</sub> + h</sup>p(t)dt
+ * </pre>
+ * where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on
+ * the derivatives at previous steps f<sub>n-k+1</sub>, f<sub>n-k+2</sub> ...
+ * f<sub>n</sub>, f<sub>n</sub> and f<sub>n+1</sub>.</p>
+ *
+ * @see AdamsMoultonIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+class AdamsMoultonStepInterpolator extends MultistepStepInterpolator {
+
+    /** Serializable version identifier */
+    private static final long serialVersionUID = 735568489801241899L;
+
+    /** Neville's interpolation array. */
+    private double[] neville;
+
+    /** Integration rollback array. */
+    private double[] rollback;
+
+    /** &gamma; star array. */
+    private double[] gammaStar;
+
+    /** Backward differences array. */
+    private int[][] bdArray;
+
+    /** Original non-truncated step end time. */
+    private double nonTruncatedEnd;
+
+    /** Original non-truncated step size. */
+    private double nonTruncatedH;
+
+    /** Simple constructor.
+     * This constructor builds an instance that is not usable yet, the
+     * {@link AbstractStepInterpolator#reinitialize} method should be called
+     * before using the instance in order to initialize the internal arrays. This
+     * constructor is used only in order to delay the initialization in
+     * some cases.
+     */
+    public AdamsMoultonStepInterpolator() {
+    }
+
+    /** Copy constructor.
+     * @param interpolator interpolator to copy from. The copy is a deep
+     * copy: its arrays are separated from the original arrays of the
+     * instance
+     */
+    public AdamsMoultonStepInterpolator(final AdamsMoultonStepInterpolator interpolator) {
+        super(interpolator);
+        nonTruncatedEnd = interpolator.nonTruncatedEnd;
+        nonTruncatedH   = interpolator.nonTruncatedH;
+    }
+
+    /** {@inheritDoc} */
+    protected StepInterpolator doCopy() {
+        return new AdamsMoultonStepInterpolator(this);
+    }
+
+    /** {@inheritDoc} */
+    protected void initializeCoefficients() {
+
+        neville  = new double[previousF.length];
+        rollback = new double[previousF.length];
+
+        bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
+
+        Fraction[] fGammaStar = AdamsMoultonIntegrator.computeGammaStarArray(previousF.length);
+        gammaStar = new double[fGammaStar.length];
+        for (int i = 0; i < fGammaStar.length; ++i) {
+            gammaStar[i] = fGammaStar[i].doubleValue();
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public void storeTime(final double t) {
+        nonTruncatedEnd = t;
+        nonTruncatedH   = nonTruncatedEnd - previousTime;
+        super.storeTime(t);
+    }
+
+    /** Truncate a step.
+     * <p>Truncating a step is necessary when an event is triggered
+     * before the nominal end of the step.</p>
+     */
+    void truncateStep(final double truncatedEndTime) {
+        currentTime = truncatedEndTime;
+        h = currentTime - previousTime;
+    }
+
+    /** {@inheritDoc} */
+    public void setInterpolatedTime(final double time)
+        throws DerivativeException {
+        interpolatedTime = time;
+        final double oneMinusThetaH = nonTruncatedEnd - interpolatedTime;
+        final double theta = (nonTruncatedH == 0) ?
+                             0 : (nonTruncatedH - oneMinusThetaH) / nonTruncatedH;
+        computeInterpolatedState(theta, oneMinusThetaH);
+    }
+
+    /** {@inheritDoc} */
+    protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
+        interpolateDerivatives();
+        interpolateState(theta);
+    }
+
+    /** Interpolate the derivatives.
+     * <p>The Adams method is based on a polynomial interpolation of the
+     * derivatives based on the preceding steps. So the interpolation of
+     * the derivatives here is strictly equivalent: it is a simple polynomial
+     * interpolation.</p>
+     */
+    private void interpolateDerivatives() {
+
+        for (int i = 0; i < interpolatedDerivatives.length; ++i) {
+
+            // initialize the Neville's interpolation algorithm
+            for (int k = 0; k < previousF.length; ++k) {
+                neville[k] = previousF[k][i];
+            }
+
+            // combine the contributions of each points
+            for (int l = 1; l < neville.length; ++l) {
+                for (int m = neville.length - 1; m >= l; --m) {
+                    final double xm   = previousT[m];
+                    final double xmMl = previousT[m - l];
+                    neville[m] = ((interpolatedTime - xm) * neville[m-1] +
+                                  (xmMl - interpolatedTime) * neville[m]) / (xmMl - xm);
+                }
+            }
+
+            // the interpolation polynomial value is in the array last element
+            interpolatedDerivatives[i] = neville[neville.length - 1];
+
+        }
+
+    }
+
+    /** Interpolate the state.
+     * <p>The Adams method is based on a polynomial interpolation of the
+     * derivatives based on the preceding steps. The polynomial model is
+     * integrated analytically throughout the last step. Using the notations
+     * found in the second edition of the first volume (Nonstiff Problems)
+     * of the reference book by Hairer, Norsett and Wanner: <i>Solving Ordinary
+     * Differential Equations</i> (Springer-Verlag, ISBN 3-540-56670-8), this
+     * process leads to the following expression:</p>
+     * <pre>
+     * y<sub>n+1</sub> = y<sub>n</sub> +
+     * h &times; &sum;<sub>j=0</sub><sup>j=k</sup> &gamma;<sub>j</sub><sup>*</sup>&nabla;<sup>j</sup>f<sub>n+1</sub>
+     * </pre>
+     * <p>In the previous expression, the &gamma;<sub>j</sub><sup>*</sup> terms are the
+     * ones that result from the analytical integration, and can be computed form
+     * the binomial coefficients C<sub>j</sub><sup>-s</sup>:</p>
+     * <p>
+     * &gamma;<sub>j</sub><sup>*</sup> = (-1)<sup>j</sup>&int;<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>1-s</sup>ds
+     * </p>
+     * <p>In order to interpolate the state in a manner that is consistent with the
+     * integration scheme, we simply subtract from the current state (at the end of the step)
+     * the integral computed from interpolation time to step end time.</p>
+     * <p>
+     * &eta;<sub>j</sub><sup>*</sup>(&theta;)=
+     * (-1)<sup>j</sup>&int;<sub>&theta;</sub><sup>1</sup>C<sub>j</sub><sup>1-s</sup>ds
+     * </p>
+     * The method described in the Hairer, Norsett and Wanner book to compute &gamma;<sub>j</sub><sup>*</sup>
+     * is easily extended to compute &gamma;<sub>j</sub><sup>*</sup>(&theta;)=
+     * (-1)<sup>j</sup>&int;<sub>0</sub><sup>&theta;</sup>C<sub>j</sub><sup>1-s</sup>ds. From this,
+     * we can compute &eta;<sub>j</sub><sup>*</sup>(&theta;) =
+     * &gamma;<sub>j</sub><sup>*</sup>-&gamma;<sub>j</sub><sup>*</sup>(&theta;).
+     * The first few values are:</p>
+     * <table>
+     * <tr><td>j</td><td>&gamma;<sub>j</sub><sup>*</sup></td><td>&gamma;<sub>j</sub><sup>*</sup>(&theta;)</td><td>&eta;<sub>j</sub><sup>*</sup>(&theta;)</td></tr>
+     * <tr><td>0</td><td>1</td><td>&theta;</td><td>1-&theta;</td></tr>
+     * <tr><td>1</td><td>-1/2</td><td>(&theta;<sup>2</sup>-2&theta;)/2</td><td>(-1+2&theta;-&theta;<sup>2</sup>)/2</td></tr>
+     * <tr><td>2</td><td>-1/12</td><td>(2&theta;<sup>3</sup>-3&theta;<sup>2</sup>)/12</td><td>(-1+3&theta;<sup>2</sup>-2&theta;<sup>3</sup>)/12</td></tr>
+     * </table>
+     * <p>
+     * The &eta;<sub>j</sub>(&theta;) functions appear to be polynomial ones. As expected,
+     * we see that &eta;<sub>j</sub>(1)= 0. The recurrence relation derived for
+     * &gamma;<sub>j</sub>(&theta;) is:
+     * </p>
+     * <p>
+     * &sum<sub>j=0</sub><sup>j=m</sup>&gamma;<sub>j</sub><sup>*</sup>(&theta;)/(m+1-j) =
+     * 1/(m+1)! &prod;<sub>k=0</sub><sup>k=m</sup>(&theta;+k-1)
+     * </p>
+     * @param theta location of the interpolation point within the last step
+     */
+    private void interpolateState(final double theta) {
+
+        // compute the integrals to remove from the final state
+        computeRollback(previousT.length - 1, theta);
+
+        // remove these integrals from the final state
+        for (int j = 0; j < interpolatedState.length; ++j) {
+            double sum = 0;
+            for (int l = 0; l < previousT.length; ++l) {
+                sum += rollback[l] * previousF[l][j];
+            }
+            interpolatedState[j] = currentState[j] - h * sum;
+        }
+
+    }
+
+    /** Compute the rollback coefficients.
+     * @param order order of the integration method
+     * @param theta current value for theta
+     */
+    private void computeRollback(final int order, final double theta) {
+
+        // compute the gamma star(theta) values from the recurrence relation
+        double product = theta - 1;
+        rollback[0]  = theta;
+        for (int i = 1; i <= order; ++i) {
+            product *= (i - 1 + theta) / (i + 1);
+            double gStar = product;
+            for (int j = 1; j <= i; ++j) {
+                gStar -= rollback[i - j] / (j + 1);
+            }
+            rollback[i] = gStar;
+        }
+
+        // subtract it from gamma star to get eta star(theta)
+        for (int i = 0; i <= order; ++i) {
+            rollback[i] -= gammaStar[i];
+        }
+
+        // combine the eta star integrals with the backward differences array
+        // to get the rollback coefficients
+        for (int i = 0; i <= order; ++i) {
+            double f = 0;
+            for (int j = i; j <= order; ++j) {
+                f -= rollback[j] * bdArray[j][i];
+            }
+            rollback[i] = f;
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public void writeExternal(final ObjectOutput out)
+        throws IOException {
+        super.writeExternal(out);
+        out.writeDouble(nonTruncatedEnd);
+    }
+
+    /** {@inheritDoc} */
+    public void readExternal(final ObjectInput in)
+        throws IOException {
+        nonTruncatedEnd = in.readDouble();
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,325 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import java.util.Arrays;
+
+import org.apache.commons.math.ode.AbstractIntegrator;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.FirstOrderIntegrator;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.ODEIntegrator;
+import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.events.EventException;
+import org.apache.commons.math.ode.events.EventHandler;
+import org.apache.commons.math.ode.events.EventState;
+import org.apache.commons.math.ode.sampling.FixedStepHandler;
+import org.apache.commons.math.ode.sampling.StepHandler;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+import org.apache.commons.math.ode.sampling.StepNormalizer;
+
+/**
+ * This class is the base class for multistep integrators for Ordinary
+ * Differential Equations.
+ *
+ * @see AdamsBashforthIntegrator
+ * @see AdamsMoultonIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public abstract class MultistepIntegrator extends AbstractIntegrator {
+
+    /** Starter integrator. */
+    private FirstOrderIntegrator starter;
+
+    /** Previous steps times. */
+    protected double[] previousT;
+
+    /** Previous steps derivatives. */
+    protected double[][] previousF;
+
+    /** Time of last detected reset. */
+    private double resetTime;
+
+    /** Prototype of the step interpolator. */
+    protected MultistepStepInterpolator prototype;
+                                           
+    /**
+     * Build a multistep integrator with the given number of steps.
+     * <p>The default starter integrator is set to the {@link
+     * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
+     * some defaults settings.</p>
+     * @param name name of the method
+     * @param k number of steps of the multistep method
+     * (including the one being computed)
+     * @param prototype prototype of the step interpolator to use
+     */
+    protected MultistepIntegrator(final String name, final int k,
+                                  final MultistepStepInterpolator prototype) {
+        super(name);
+        starter = new DormandPrince853Integrator(1.0e-6, 1.0e6, 1.0e-5, 1.0e-6);
+        previousT = new double[k];
+        previousF = new double[k][];
+        this.prototype = prototype;
+    }
+
+    /**
+     * Get the starter integrator.
+     * @return starter integrator
+     */
+    public ODEIntegrator getStarterIntegrator() {
+        return starter;
+    }
+
+    /**
+     * Set the starter integrator.
+     * <p>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.</p>
+     * @param starter starter integrator
+     */
+    public void setStarterIntegrator(FirstOrderIntegrator starter) {
+        this.starter = starter;
+    }
+
+    /** Start the integration.
+     * <p>This method computes the first few steps of the multistep method,
+     * using the underlying starter integrator, ensuring the returned steps
+     * all belong to the same smooth range.</p>
+     * <p>In order to ensure smoothness, the start phase is automatically
+     * restarted when a state or derivative reset is triggered by the
+     * registered events handlers before this start phase is completed. As
+     * an example, consider integrating a differential equation from t=0
+     * to t=100 with a 4 steps method and step size equal to 0.2. If an event
+     * resets the state at t=0.5, the start phase will not end at t=0.7 with
+     * steps at [0.0, 0.2, 0.4, 0.6] but instead will end at t=1.1 with steps
+     * at [0.5, 0.7, 0.9, 1.1].</p>
+     * <p>A side effect of the need for smoothness is that an ODE triggering
+     * short period regular resets will remain in the start phase throughout
+     * the integration range if the step size or the number of steps to store
+     * are too large.</p>
+     * <p>If the start phase ends prematurely (because of some triggered event
+     * for example), then the time of latest previous steps will be set to
+     * <code>Double.NaN</code>.</p>
+     * @param n number of steps to store
+     * @param h signed step size to use for the first steps
+     * @param manager discrete events manager to use
+     * @param equations differential equations to integrate
+     * @param t0 initial time
+     * @param y state vector: contains the initial value of the state vector at t0,
+     * will be used to put the state vector at each successful step and hence
+     * contains the final value at the end of the start phase
+     * @return time of the end of the start phase
+     * @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 double start(final int n, final double h,
+                           final CombinedEventsManager manager,
+                           final FirstOrderDifferentialEquations equations,
+                           final double t0, final double[] y)
+        throws DerivativeException, IntegratorException {
+
+        // clear the first steps
+        Arrays.fill(previousT, Double.NaN);
+        Arrays.fill(previousF, null);
+
+        // configure the event handlers
+        starter.clearEventHandlers();
+        for (EventState state : manager.getEventsStates()) {
+            starter.addEventHandler(new ResetCheckingWrapper(state.getEventHandler()),
+                                    state.getMaxCheckInterval(),
+                                    state.getConvergence(), state.getMaxIterationCount());
+        }
+
+        // configure the step handlers
+        starter.clearStepHandlers();
+        for (final StepHandler handler : stepHandlers) {
+            // add the user defined step handlers, filtering out the isLast indicator
+            starter.addStepHandler(new FilteringWrapper(handler));
+        }
+
+        // add one specific step handler to store the first steps
+        final StoringStepHandler store = new StoringStepHandler(n);
+        starter.addStepHandler(new StepNormalizer(h, store));
+
+        // integrate over the first few steps, ensuring no intermediate reset occurs
+        double t = t0;
+        double stopTime = Double.NaN;
+        do {
+            resetTime = Double.NaN;
+            store.restart();
+            // we overshoot by 1/10000 step the end to make sure we get don't miss the last point
+            stopTime = starter.integrate(equations, t, y, t + (n - 0.9999) * h, y);
+            if (!Double.isNaN(resetTime)) {
+                // there was an intermediate reset, we restart
+                t = resetTime;
+            }
+        } while (!Double.isNaN(resetTime));
+
+        // clear configuration
+        starter.clearEventHandlers();
+        starter.clearStepHandlers();
+
+        if (store.getFinalState() != null) {
+            System.arraycopy(store.getFinalState(), 0, y, 0, y.length);
+        }
+        return stopTime;
+
+    }
+
+    /** Rotate the previous steps arrays.
+     */
+    protected void rotatePreviousSteps() {
+        final double[] rolled = previousF[previousT.length - 1];
+        for (int k = previousF.length - 1; k > 0; --k) {
+            previousT[k] = previousT[k - 1];
+            previousF[k] = previousF[k - 1];
+        }
+        previousF[0] = rolled;
+    }
+
+    /** Event handler wrapper to check if state or derivatives have been reset. */
+    private class ResetCheckingWrapper implements EventHandler {
+
+        /** Serializable version identifier. */
+        private static final long serialVersionUID = 4922660285376467937L;
+
+        /** Wrapped event handler. */
+        private final EventHandler handler;
+
+        /** Build a new instance.
+         * @param handler event handler to wrap
+         */
+        public ResetCheckingWrapper(final EventHandler handler) {
+            this.handler = handler;
+        }
+
+        /** {@inheritDoc} */
+        public int eventOccurred(double t, double[] y) throws EventException {
+            final int action = handler.eventOccurred(t, y);
+            if ((action == RESET_DERIVATIVES) || (action == RESET_STATE)) {
+                // a singularity has been encountered
+                // we need to restart the start phase
+                resetTime = t;
+                return STOP;
+            }
+            return action;
+        }
+
+        /** {@inheritDoc} */
+        public double g(double t, double[] y) throws EventException {
+            return handler.g(t, y);
+        }
+
+        /** {@inheritDoc} */
+        public void resetState(double t, double[] y) throws EventException {
+            handler.resetState(t, y);
+        }
+        
+    }
+
+    /** Step handler wrapper filtering out the isLast indicator. */
+    private class FilteringWrapper implements StepHandler {
+
+        /** Serializable version identifier. */
+        private static final long serialVersionUID = 4607975253344802232L;
+
+        /** Wrapped step handler. */
+        private final StepHandler handler;
+
+        /** Build a new instance.
+         * @param handler step handler to wrap
+         */
+        public FilteringWrapper(final StepHandler handler) {
+            this.handler = handler;
+        }
+
+        /** {@inheritDoc} */
+        public void handleStep(StepInterpolator interpolator, boolean isLast)
+                throws DerivativeException {
+            // we force the isLast indicator to false EXCEPT if some event handler triggered a stop
+            handler.handleStep(interpolator, eventsHandlersManager.stop());
+        }
+
+        /** {@inheritDoc} */
+        public boolean requiresDenseOutput() {
+            return handler.requiresDenseOutput();
+        }
+
+        /** {@inheritDoc} */
+        public void reset() {
+            handler.reset();
+        }
+        
+    }
+
+    /** Specialized step handler storing the first few steps. */
+    private class StoringStepHandler implements FixedStepHandler {
+
+        /** Serializable version identifier. */
+        private static final long serialVersionUID = 4592974435520688797L;
+
+        /** Number of steps to store. */
+        private final int n;
+
+        /** Counter for already stored steps. */
+        private int count;
+
+        /** Final state. */
+        private double[] finalState;
+
+        /** Build a new instance.
+         * @param number of steps to store
+         */
+        public StoringStepHandler(final int n) {
+            this.n = n;
+            restart();
+        }
+
+        /** Restart storage.
+         */
+        public void restart() {
+            count = 0;
+            finalState = null;
+        }
+
+        /** Get the final state.
+         * @return final state
+         */
+        public double[] getFinalState() {
+            return finalState;
+        }
+
+        /** {@inheritDoc} */
+        public void handleStep(final double t, final double[] y, final double[] yDot,
+                               final boolean isLast) {
+            if (count++ < n) {
+                previousT[n - count] = t;
+                previousF[n - count] = yDot.clone();
+                if (count == n) {
+                    finalState = y.clone();
+                }
+            }
+        }
+
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,165 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+
+/** This class represents an interpolator over the last step during an
+ * ODE integration for multistep integrators.
+ *
+ * @see MultistepIntegrator
+ *
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+abstract class MultistepStepInterpolator
+    extends AbstractStepInterpolator {
+
+    /** Previous steps times. */
+    protected double[] previousT;
+
+    /** Previous steps derivatives. */
+    protected double[][] previousF;
+
+    /** Simple constructor.
+     * This constructor builds an instance that is not usable yet, the
+     * {@link #reinitialize} method should be called before using the
+     * instance in order to initialize the internal arrays. This
+     * constructor is used only in order to delay the initialization in
+     * some cases. The {@link MultistepIntegrator} classe uses the
+     * prototyping design pattern to create the step interpolators by
+     * cloning an uninitialized model and latter initializing the copy.
+     */
+    protected MultistepStepInterpolator() {
+        previousT = null;
+        previousF = null;
+    }
+
+    /** Copy constructor.
+
+     * <p>The copied interpolator should have been finalized before the
+     * copy, otherwise the copy will not be able to perform correctly any
+     * interpolation and will throw a {@link NullPointerException}
+     * later. Since we don't want this constructor to throw the
+     * exceptions finalization may involve and since we don't want this
+     * method to modify the state of the copied interpolator,
+     * finalization is <strong>not</strong> done automatically, it
+     * remains under user control.</p>
+
+     * <p>The copy is a deep copy: its arrays are separated from the
+     * original arrays of the instance.</p>
+
+     * @param interpolator interpolator to copy from.
+
+     */
+    public MultistepStepInterpolator(final MultistepStepInterpolator interpolator) {
+
+        super(interpolator);
+
+        if (interpolator.currentState != null) {
+            previousT = interpolator.previousT.clone();
+            previousF = new double[interpolator.previousF.length][];
+            for (int k = 0; k < interpolator.previousF.length; ++k) {
+                previousF[k] = interpolator.previousF[k].clone();
+            }
+            initializeCoefficients();
+        } else {
+            previousT = null;
+            previousF = null;
+        }
+
+    }
+
+    /** Reinitialize the instance
+     * @param y reference to the integrator array holding the state at
+     * the end of the step
+     * @param previousT reference to the integrator array holding the times
+     * of the previous steps
+     * @param previousF reference to the integrator array holding the
+     * previous slopes
+     * @param forward integration direction indicator
+     */
+    public void reinitialize(final double[] y,
+                             final double[] previousT, final double[][] previousF,
+                             final boolean forward) {
+        reinitialize(y, forward);
+        this.previousT = previousT;
+        this.previousF = previousF;
+        initializeCoefficients();
+    }
+
+    /** Initialize the coefficients arrays.
+     */
+    protected abstract void initializeCoefficients();
+
+    /** {@inheritDoc} */
+    public void writeExternal(final ObjectOutput out)
+    throws IOException {
+
+        // save the state of the base class
+        writeBaseExternal(out);
+
+        // save the local attributes
+        out.writeInt(previousT.length);
+        for (int k = 0; k < previousF.length; ++k) {
+            out.writeDouble(previousT[k]);
+            for (int i = 0; i < currentState.length; ++i) {
+                out.writeDouble(previousF[k][i]);
+            }
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public void readExternal(final ObjectInput in)
+    throws IOException {
+
+        // read the base class 
+        final double t = readBaseExternal(in);
+
+        // read the local attributes
+        final int kMax = in.readInt();
+        previousT = new double[kMax];
+        previousF = new double[kMax][];
+        for (int k = 0; k < kMax; ++k) {
+            previousT[k] = in.readDouble();
+            previousF[k] = new double[currentState.length];
+            for (int i = 0; i < currentState.length; ++i) {
+                previousF[k][i] = in.readDouble();
+            }
+        }
+
+        // initialize the coefficients
+        initializeCoefficients();
+
+        try {
+            // we can now set the interpolated time and state
+            setInterpolatedTime(t);
+        } catch (DerivativeException e) {
+            throw new IOException(e.getMessage());
+        }
+
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml?rev=676737&r1=676736&r2=676737&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml Mon Jul 14 14:13:23 2008
@@ -39,6 +39,13 @@
   </properties>
   <body>
     <release version="2.0" date="TBD" description="TBD">
+      <action dev="luc" type="add">
+        New ODE integrators have been added: the explicit Adams-Bashforth and implicit
+        Adams-Moulton multistep methods. These methods support customizable starter
+        integrators and support discrete events even during the start phase.
+        All these methods provide the same rich features has the existing ones:
+        continuous output, step handlers, discrete events, G-stop ...
+      </action>
       <action dev="luc" type="fix" issue="MATH-214" >
         Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by
         a truncation of the last step only.

Added: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,206 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.nonstiff;
+
+import junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderIntegrator;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.events.EventHandler;
+
+public class AdamsBashforthIntegratorTest
+  extends TestCase {
+
+  public AdamsBashforthIntegratorTest(String name) {
+    super(name);
+  }
+
+  public void testCoefficients() {
+
+      double[] coeffs1 = new AdamsBashforthIntegrator(1, 0.01).getCoeffs();
+      assertEquals(1, coeffs1.length);
+      assertEquals(1.0, coeffs1[0], 1.0e-16);
+
+      double[] coeffs2 = new AdamsBashforthIntegrator(2, 0.01).getCoeffs();
+      assertEquals(2, coeffs2.length);
+      assertEquals( 3.0 / 2.0, coeffs2[0], 1.0e-16);
+      assertEquals(-1.0 / 2.0, coeffs2[1], 1.0e-16);
+
+      double[] coeffs3 = new AdamsBashforthIntegrator(3, 0.01).getCoeffs();
+      assertEquals(3, coeffs3.length);
+      assertEquals( 23.0 / 12.0, coeffs3[0], 1.0e-16);
+      assertEquals(-16.0 / 12.0, coeffs3[1], 1.0e-16);
+      assertEquals(  5.0 / 12.0, coeffs3[2], 1.0e-16);
+
+      double[] coeffs4 = new AdamsBashforthIntegrator(4, 0.01).getCoeffs();
+      assertEquals(4, coeffs4.length);
+      assertEquals( 55.0 / 24.0, coeffs4[0], 1.0e-16);
+      assertEquals(-59.0 / 24.0, coeffs4[1], 1.0e-16);
+      assertEquals( 37.0 / 24.0, coeffs4[2], 1.0e-16);
+      assertEquals( -9.0 / 24.0, coeffs4[3], 1.0e-16);
+
+      double[] coeffs5 = new AdamsBashforthIntegrator(5, 0.01).getCoeffs();
+      assertEquals(5, coeffs5.length);
+      assertEquals( 1901.0 / 720.0, coeffs5[0], 1.0e-16);
+      assertEquals(-2774.0 / 720.0, coeffs5[1], 1.0e-16);
+      assertEquals( 2616.0 / 720.0, coeffs5[2], 1.0e-16);
+      assertEquals(-1274.0 / 720.0, coeffs5[3], 1.0e-16);
+      assertEquals(  251.0 / 720.0, coeffs5[4], 1.0e-16);
+
+      double[] coeffs6 = new AdamsBashforthIntegrator(6, 0.01).getCoeffs();
+      assertEquals(6, coeffs6.length);
+      assertEquals( 4277.0 / 1440.0, coeffs6[0], 1.0e-16);
+      assertEquals(-7923.0 / 1440.0, coeffs6[1], 1.0e-16);
+      assertEquals( 9982.0 / 1440.0, coeffs6[2], 1.0e-16);
+      assertEquals(-7298.0 / 1440.0, coeffs6[3], 1.0e-16);
+      assertEquals( 2877.0 / 1440.0, coeffs6[4], 1.0e-16);
+      assertEquals( -475.0 / 1440.0, coeffs6[5], 1.0e-16);
+
+      double[] coeffs7 = new AdamsBashforthIntegrator(7, 0.01).getCoeffs();
+      assertEquals(7, coeffs7.length);
+      assertEquals( 198721.0 / 60480.0, coeffs7[0], 1.0e-16);
+      assertEquals(-447288.0 / 60480.0, coeffs7[1], 1.0e-16);
+      assertEquals( 705549.0 / 60480.0, coeffs7[2], 1.0e-16);
+      assertEquals(-688256.0 / 60480.0, coeffs7[3], 1.0e-16);
+      assertEquals( 407139.0 / 60480.0, coeffs7[4], 1.0e-16);
+      assertEquals(-134472.0 / 60480.0, coeffs7[5], 1.0e-16);
+      assertEquals(  19087.0 / 60480.0, coeffs7[6], 1.0e-16);
+
+      double[] coeffs8 = new AdamsBashforthIntegrator(8, 0.01).getCoeffs();
+      assertEquals(8, coeffs8.length);
+      assertEquals(  434241.0 / 120960.0, coeffs8[0], 1.0e-16);
+      assertEquals(-1152169.0 / 120960.0, coeffs8[1], 1.0e-16);
+      assertEquals( 2183877.0 / 120960.0, coeffs8[2], 1.0e-16);
+      assertEquals(-2664477.0 / 120960.0, coeffs8[3], 1.0e-16);
+      assertEquals( 2102243.0 / 120960.0, coeffs8[4], 1.0e-16);
+      assertEquals(-1041723.0 / 120960.0, coeffs8[5], 1.0e-16);
+      assertEquals(  295767.0 / 120960.0, coeffs8[6], 1.0e-16);
+      assertEquals(  -36799.0 / 120960.0, coeffs8[7], 1.0e-16);
+
+      double[] coeffs9 = new AdamsBashforthIntegrator(9, 0.01).getCoeffs();
+      assertEquals(9, coeffs9.length);
+      assertEquals(  14097247.0 / 3628800.0, coeffs9[0], 1.0e-16);
+      assertEquals( -43125206.0 / 3628800.0, coeffs9[1], 1.0e-16);
+      assertEquals(  95476786.0 / 3628800.0, coeffs9[2], 1.0e-16);
+      assertEquals(-139855262.0 / 3628800.0, coeffs9[3], 1.0e-16);
+      assertEquals( 137968480.0 / 3628800.0, coeffs9[4], 1.0e-16);
+      assertEquals( -91172642.0 / 3628800.0, coeffs9[5], 1.0e-16);
+      assertEquals(  38833486.0 / 3628800.0, coeffs9[6], 1.0e-16);
+      assertEquals(  -9664106.0 / 3628800.0, coeffs9[7], 1.0e-16);
+      assertEquals(   1070017.0 / 3628800.0, coeffs9[8], 1.0e-16);
+
+  }
+
+  public void testDimensionCheck() {
+    try  {
+      TestProblem1 pb = new TestProblem1();
+      new AdamsBashforthIntegrator(3, 0.01).integrate(pb,
+                                                      0.0, new double[pb.getDimension()+10],
+                                                      1.0, new double[pb.getDimension()+10]);
+        fail("an exception should have been thrown");
+    } catch(DerivativeException de) {
+      fail("wrong exception caught");
+    } catch(IntegratorException ie) {
+    }
+  }
+  
+  public void testDecreasingSteps()
+    throws DerivativeException, IntegratorException {
+
+    TestProblemAbstract[] problems = TestProblemFactory.getProblems();
+    for (int k = 0; k < problems.length; ++k) {
+
+      double previousError = Double.NaN;
+      for (int i = 6; i < 10; ++i) {
+
+        TestProblemAbstract pb  = (TestProblemAbstract) problems[k].clone();
+        double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i);
+
+        FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, step);
+        TestProblemHandler handler = new TestProblemHandler(pb, integ);
+        integ.addStepHandler(handler);
+        EventHandler[] functions = pb.getEventsHandlers();
+        for (int l = 0; l < functions.length; ++l) {
+          integ.addEventHandler(functions[l],
+                                Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
+        }
+        double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+                                          pb.getFinalTime(), new double[pb.getDimension()]);
+        if (functions.length == 0) {
+          assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
+        }
+
+        double error = handler.getMaximalValueError();
+        if (i > 6) {
+          assertTrue(error < Math.abs(previousError));
+        }
+        previousError = error;
+
+      }
+
+    }
+
+  }
+
+  public void testSmallStep()
+    throws DerivativeException, IntegratorException {
+
+    TestProblem1 pb  = new TestProblem1();
+    double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+
+    FirstOrderIntegrator integ = new AdamsBashforthIntegrator(3, step);
+    TestProblemHandler handler = new TestProblemHandler(pb, integ);
+    integ.addStepHandler(handler);
+    integ.integrate(pb,
+                    pb.getInitialTime(), pb.getInitialState(),
+                    pb.getFinalTime(), new double[pb.getDimension()]);
+
+   assertTrue(handler.getLastError() < 2.0e-9);
+   assertTrue(handler.getMaximalValueError() < 3.0e-8);
+   assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+   assertEquals("Adams-Bashforth", integ.getName());
+
+  }
+
+  public void testBigStep()
+    throws DerivativeException, IntegratorException {
+
+    TestProblem1 pb  = new TestProblem1();
+    double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
+
+    FirstOrderIntegrator integ = new AdamsBashforthIntegrator(3, step);
+    TestProblemHandler handler = new TestProblemHandler(pb, integ);
+    integ.addStepHandler(handler);
+    integ.integrate(pb,
+                    pb.getInitialTime(), pb.getInitialState(),
+                    pb.getFinalTime(), new double[pb.getDimension()]);
+
+    assertTrue(handler.getLastError() > 0.05);
+    assertTrue(handler.getMaximalValueError() > 0.1);
+    assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+
+  }
+  
+  public static Test suite() {
+    return new TestSuite(AdamsBashforthIntegratorTest.class);
+  }
+
+}

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision



Mime
View raw message