commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r786877 [1/2] - in /commons/proper/math/trunk/src: java/org/apache/commons/math/ode/ java/org/apache/commons/math/ode/nonstiff/ java/org/apache/commons/math/ode/sampling/ test/org/apache/commons/math/ode/nonstiff/ test/org/apache/commons/ma...
Date Sat, 20 Jun 2009 18:15:54 GMT
Author: luc
Date: Sat Jun 20 18:15:54 2009
New Revision: 786877

URL: http://svn.apache.org/viewvc?rev=786877&view=rev
Log:
changed the Adams-Bashforth and Adams-Moulton multistep integrators to adaptive stepsize.
this was made possible thanks to the Nordsieck representation of integration state

Added:
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java   (with props)
Modified:
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
    commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
    commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
    commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
    commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java

Added: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java?rev=786877&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java Sat Jun 20 18:15:54 2009
@@ -0,0 +1,326 @@
+/*
+ * 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;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.apache.commons.math.fraction.BigFraction;
+import org.apache.commons.math.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math.linear.Array2DRowRealMatrix;
+import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
+import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
+import org.apache.commons.math.linear.FieldDecompositionSolver;
+import org.apache.commons.math.linear.FieldLUDecompositionImpl;
+import org.apache.commons.math.linear.FieldMatrix;
+import org.apache.commons.math.linear.MatrixUtils;
+import org.apache.commons.math.linear.MatrixVisitorException;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator;
+import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
+
+/** Transformer for Nordsieck vectors.
+ * <p>This class i used by {@link MultistepIntegrator multistep integrators}
+ * to convert between classical representation with several previous first
+ * derivatives and Nordsieck representation with higher order scaled derivatives.</p>
+ *
+ * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
+ * <pre>
+ * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
+ * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
+ * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
+ * ...
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
+ * </pre></p>
+ *
+ * <p>With the previous definition, the classical representation of multistep methods
+ * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
+ * q<sub>n</sub> where q<sub>n</sub> is defined as:
+ * <pre>
+ *   q<sub>n</sub> = [ s<sub>1</sub>(n-1) s<sub>1</sub>(n-2) ... s<sub>1</sub>(n-(k-1)) ]<sup>T</sup>
+ * </pre>
+ * (we omit the k index in the notation for clarity).</p>
+ *
+ * <p>Another possible representation uses the Nordsieck vector with
+ * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
+ * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
+ * </pre>
+ * (here again we omit the k index in the notation for clarity)
+ * </p>
+ *
+ * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
+ * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
+ * for degree k polynomials.
+ * <pre>
+ * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + &sum;<sub>j</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
+ * </pre>
+ * The previous formula can be used with several values for i to compute the transform between
+ * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
+ * and q<sub>n</sub> resulting from the Taylor series formulas above is:
+ * <pre>
+ * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
+ * </pre>
+ * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)&times;(k-1) matrix built
+ * with the j (-i)<sup>j-1</sup> terms:
+ * <pre>
+ *        [  -2   3   -4    5  ... ]
+ *        [  -4  12  -32   80  ... ]
+ *   P =  [  -6  27 -108  405  ... ]
+ *        [  -8  48 -256 1280  ... ]
+ *        [          ...           ]
+ * </pre></p>
+ *
+ * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
+ * classical representation and Nordsieck vector at step start. The resulting Q matrix is simply
+ * the absolute value of matrix P.</p>
+ * 
+ * <p>Using the Nordsieck vector has several advantages:
+ * <ul>
+ *   <li>it greatly simplifies step interpolation as the interpolator mainly applies
+ *   Taylor series formulas,</li>
+ *   <li>it simplifies step changes that occur when discrete events that truncate
+ *   the step are triggered,</li>
+ *   <li>it allows to extend the methods in order to support adaptive stepsize (not implemented yet).</li>
+ * </ul></p>
+ * 
+ * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ *   <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ *   <li>r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * where A is a rows shifting matrix (the lower left part is an identity matrix):
+ * <pre>
+ *        [ 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 ]
+ * </pre></p>
+ *
+ * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ *   <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ *   <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
+ *   <li>R<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * where A is a rows shifting matrix (the lower left part is an identity matrix):
+ * <pre>
+ *        [ 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 ]
+ * </pre>
+ * From this predicted vector, the corrected vector is computed as follows:
+ * <ul>
+ *   <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... &plusmn;1 ] r<sub>n+1</sub></li>
+ *   <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ *   <li>r<sub>n+1</sub> = R<sub>n+1</sub> + (s<sub>1</sub>(n+1) - S<sub>1</sub>(n+1)) P<sup>-1</sup> u</li>
+ * </ul>
+ * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
+ * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
+ * represent the corrected states.</p>
+ *
+ * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
+ * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
+ * they only depend on k. This class handles these transformations.</p>
+ *
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class NordsieckTransformer {
+
+    /** Cache for already computed coefficients. */
+    private static final Map<Integer, NordsieckTransformer> cache =
+        new HashMap<Integer, NordsieckTransformer>();
+
+    /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
+    private final RealMatrix initialization;
+
+    /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
+    private final RealMatrix update;
+
+    /** Update coefficients of the higher order derivatives wrt y'. */
+    private final double[] c1;
+
+    /** Simple constructor.
+     * @param nSteps number of steps of the multistep method
+     * (including the one being computed)
+     */
+    private NordsieckTransformer(final int nSteps) {
+
+        // compute exact coefficients
+        FieldMatrix<BigFraction> bigP = buildP(nSteps);
+        FieldDecompositionSolver<BigFraction> pSolver =
+            new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
+        BigFraction[] u = new BigFraction[nSteps - 1];
+        Arrays.fill(u, BigFraction.ONE);
+        BigFraction[] bigC1 = pSolver.solve(u);
+
+        // update coefficients are computed by combining transform from
+        // Nordsieck to multistep, then shifting rows to represent step advance
+        // then applying inverse transform
+        BigFraction[][] shiftedP = bigP.getData();
+        for (int i = shiftedP.length - 1; i > 0; --i) {
+            // shift rows
+            shiftedP[i] = shiftedP[i - 1];
+        }
+        shiftedP[0] = new BigFraction[nSteps - 1];
+        Arrays.fill(shiftedP[0], BigFraction.ZERO);
+        FieldMatrix<BigFraction> bigMSupdate =
+            pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
+
+        // initialization coefficients, computed from a Q matrix = abs(P)
+        bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
+            /** {@inheritDoc} */
+            @Override
+            public BigFraction visit(int row, int column, BigFraction value) {
+                return ((column & 0x1) == 0x1) ? value : value.negate();
+            }
+        });
+        FieldMatrix<BigFraction> bigQInverse =
+            new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
+
+        // convert coefficients to double
+        initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigQInverse);
+        update         = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
+        c1             = new double[nSteps - 1];
+        for (int i = 0; i < nSteps - 1; ++i) {
+            c1[i] = bigC1[i].doubleValue();
+        }
+
+    }
+
+    /** Get the Nordsieck transformer for a given number of steps.
+     * @param nSteps number of steps of the multistep method
+     * (including the one being computed)
+     * @return Nordsieck transformer for the specified number of steps
+     */
+    public static NordsieckTransformer getInstance(final int nSteps) {
+        synchronized(cache) {
+            NordsieckTransformer t = cache.get(nSteps);
+            if (t == null) {
+                t = new NordsieckTransformer(nSteps);
+                cache.put(nSteps, t);
+            }
+            return t;
+        }
+    }
+
+    /** Build the P matrix transforming multistep to Nordsieck.
+     * <p>
+     * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
+     * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
+     * The two representations share their two first components y(k) and
+     * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
+     * <pre>
+     * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
+     * </pre>
+     * </p>
+     * @param nSteps number of steps of the multistep method
+     * (including the one being computed)
+     * @return P matrix
+     */
+    private FieldMatrix<BigFraction> buildP(final int nSteps) {
+
+        final BigFraction[][] pData = new BigFraction[nSteps - 1][nSteps - 1];
+
+        for (int i = 0; i < pData.length; ++i) {
+            // build the P matrix elements from Taylor series formulas
+            final BigFraction[] pI = pData[i];
+            final int factor = -(i + 1);
+            int aj = factor;
+            for (int j = 0; j < pI.length; ++j) {
+                pI[j] = new BigFraction(aj * (j + 2));
+                aj *= factor;
+            }
+        }
+
+        return new Array2DRowFieldMatrix<BigFraction>(pData, false);
+
+    }
+
+    /** Initialize the high order scaled derivatives at step start.
+     * @param first first scaled derivative at step start
+     * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
+     * will be modified
+     * @return high order derivatives at step start
+     */
+    public RealMatrix initializeHighOrderDerivatives(final double[] first,
+                                                     final double[][] multistep) {
+        for (int i = 0; i < multistep.length; ++i) {
+            final double[] msI = multistep[i];
+            for (int j = 0; j < first.length; ++j) {
+                msI[j] -= first[j];
+            }
+        }
+        return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
+    }
+
+    /** Update the high order scaled derivatives (phase 1).
+     * <p>The complete update of high order derivatives has a form similar to:
+     * <pre>
+     * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+     * </pre>
+     * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
+     * @param highOrder high order scaled derivatives
+     * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
+     * @return updated high order derivatives
+     * @see #updateHighOrderDerivativesPhase2(double[], double[], RealMatrix)
+     */
+    public RealMatrix updateHighOrderDerivativesPhase1(final RealMatrix highOrder) {
+        return update.multiply(highOrder);
+    }
+
+    /** Update the high order scaled derivatives (phase 2).
+     * <p>The complete update of high order derivatives has a form similar to:
+     * <pre>
+     * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+     * </pre>
+     * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
+     * <p>Phase 1 of the update must already have been performed.</p>
+     * @param start first order scaled derivatives at step start
+     * @param end first order scaled derivatives at step end
+     * @param highOrder high order scaled derivatives, will be modified
+     * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
+     * @see #updateHighOrderDerivativesPhase1(RealMatrix)
+     */
+    public void updateHighOrderDerivativesPhase2(final double[] start,
+                                                 final double[] end,
+                                                 final RealMatrix highOrder) {
+        highOrder.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+            /** {@inheritDoc} */
+            @Override
+            public double visit(int row, int column, double value)
+                throws MatrixVisitorException {
+                return value + c1[row] * (start[column] - end[column]);
+            }
+        });
+    }
+
+}

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Sat Jun 20 18:15:54 2009
@@ -17,25 +17,9 @@
 
 package org.apache.commons.math.ode.nonstiff;
 
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import java.io.ObjectOutputStream;
 import java.io.Serializable;
-import java.lang.reflect.Field;
-import java.util.Arrays;
-import java.util.HashMap;
-import java.util.Map;
-
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.fraction.BigFraction;
-import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
-import org.apache.commons.math.linear.FieldLUDecompositionImpl;
-import org.apache.commons.math.linear.FieldMatrix;
-import org.apache.commons.math.linear.Array2DRowFieldMatrix;
-import org.apache.commons.math.linear.MatrixUtils;
-import org.apache.commons.math.linear.MatrixVisitorException;
+
 import org.apache.commons.math.linear.RealMatrix;
-import org.apache.commons.math.linear.Array2DRowRealMatrix;
 import org.apache.commons.math.ode.DerivativeException;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
@@ -161,48 +145,50 @@
     /** Serializable version identifier. */
     private static final long serialVersionUID = 67792782787082199L;
 
-    /** Cache for already computed coefficients. */
-    private static final Map<Integer, CachedCoefficients> cache =
-        new HashMap<Integer, CachedCoefficients>();
-
-    /** Coefficients of the method. */
-    private final transient CachedCoefficients coefficients;
-
-    /** Integration step. */
-    private final double step;
-
     /**
      * Build an Adams-Bashforth with the given order and step size.
      * @param order order of the method (must be greater than 1: due to
      * an implementation limitation the order 1 method is not supported)
-     * @param step integration step size
+     * @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 order, final double step)
+    public AdamsBashforthIntegrator(final int order,
+                                    final double minStep, final double maxStep,
+                                    final double scalAbsoluteTolerance,
+                                    final double scalRelativeTolerance)
         throws IllegalArgumentException {
+        super("Adams-Bashforth", order, order, minStep, maxStep,
+              scalAbsoluteTolerance, scalRelativeTolerance);
+    }
 
-        super("Adams-Bashforth", order);
-        if (order <= 1) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  "{0} is supported only for orders 2 or more",
-                  getName());
-        }
-
-        // cache the coefficients for each order, to avoid recomputing them
-        synchronized(cache) {
-            CachedCoefficients coeff = cache.get(order);
-            if (coeff == null) {
-                coeff = new CachedCoefficients(order);
-                cache.put(order, coeff);
-            }
-            coefficients = coeff;
-        }
-
-        this.step = Math.abs(step);
-
+    /**
+     * Build an Adams-Bashforth with the given order and step size.
+     * @param order order of the method (must be greater than 1: due to
+     * an implementation limitation the order 1 method is not supported)
+     * @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 order,
+                                    final double minStep, final double maxStep,
+                                    final double[] vecAbsoluteTolerance,
+                                    final double[] vecRelativeTolerance)
+        throws IllegalArgumentException {
+        super("Adams-Bashforth", order, order, minStep, maxStep,
+              vecAbsoluteTolerance, vecRelativeTolerance);
     }
 
     /** {@inheritDoc} */
+    @Override
     public double integrate(final FirstOrderDifferentialEquations equations,
                             final double t0, final double[] y0,
                             final double t, final double[] y)
@@ -218,66 +204,115 @@
         if (y != y0) {
             System.arraycopy(y0, 0, y, 0, n);
         }
+        final double[] yDot = new double[n];
+        final double[] yTmp = new double[y0.length];
 
         // set up an interpolator sharing the integrator arrays
         final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
         interpolator.reinitialize(y, forward);
+        final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
+        interpolatorTmp.reinitialize(yTmp, forward);
 
         // set up integration control objects
-        stepStart = t0;
-        stepSize  = forward ? step : -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, stepStart, y);
-        if (Double.isNaN(previousT[0])) {
-            return stopTime;
-        }
-        stepStart = previousT[0];
-
-        // convert to Nordsieck representation
-        double[]   scaled    = convertToNordsieckLow();
-        RealMatrix nordsieck = convertToNordsieckHigh(scaled);
+        // 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);
+        
         boolean lastStep = false;
         while (!lastStep) {
 
             // shift all data
             interpolator.shift();
 
-            // discrete events handling
-            interpolator.storeTime(stepStart + stepSize);
-            if (manager.evaluateStep(interpolator)) {
-                stepSize = manager.getEventTime() - stepStart;
+            double error = 0;
+            for (boolean loop = true; loop;) {
+
+                stepSize = hNew;
+
+                // evaluate error using the last term of the Taylor expansion
+                error = 0;
+                for (int i = 0; i < y0.length; ++i) {
+                    final double yScale = Math.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 = Math.sqrt(error / y0.length);
+
+                if (error <= 1.0) {
+
+                    // predict a first estimate of the state at step end
+                    final double stepEnd = stepStart + stepSize;
+                    interpolator.setInterpolatedTime(stepEnd);
+                    System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
+
+                    // evaluate the derivative
+                    computeDerivatives(stepEnd, yTmp, 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 RealMatrix nordsieckTmp =
+                        transformer.updateHighOrderDerivativesPhase1(nordsieck);
+                    transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
+
+                    // discrete events handling
+                    interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
+                    interpolatorTmp.storeTime(stepStart);
+                    interpolatorTmp.shift();
+                    interpolatorTmp.storeTime(stepEnd);
+                    if (manager.evaluateStep(interpolatorTmp)) {
+                        final double dt = manager.getEventTime() - stepStart;
+                        if (Math.abs(dt) <= Math.ulp(stepStart)) {
+                            // rejecting the step would lead to a too small next step, we accept it
+                            loop = false;
+                        } else {
+                            // reject the step to match exactly the next switch time
+                            hNew = dt;
+                            interpolator.rescale(hNew);
+                        }
+                    } else {
+                        // accept the step
+                        scaled    = predictedScaled;
+                        nordsieck = nordsieckTmp;
+                        interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+                        loop = false;
+                    }
+
+                } else {
+                    // 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);
+                }
+
             }
 
             // the step has been accepted (may have been truncated)
             final double nextStep = stepStart + stepSize;
+            System.arraycopy(yTmp, 0, y, 0, n);
             interpolator.storeTime(nextStep);
-            System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, n);
             manager.stepAccepted(nextStep, y);
             lastStep = manager.stop();
 
-            // update the Nordsieck vector
-            final double[] f0 = previousF[0];
-            previousT[0] = nextStep;
-            computeDerivatives(nextStep, y, f0);
-            nordsieck = coefficients.msUpdate.multiply(nordsieck);
-            final double[] end = new double[y0.length];
-            for (int j = 0; j < y0.length; ++j) {
-                end[j] = stepSize * f0[j];
-            }
-            nordsieck.walkInOptimizedOrder(new NordsieckUpdater(scaled, end, coefficients.c1));
-            scaled = end;
-            interpolator.reinitialize(nextStep, stepSize, scaled, nordsieck);
-
             // provide the step data to the step handler
             for (StepHandler handler : stepHandlers) {
+                interpolator.setInterpolatedTime(nextStep);
                 handler.handleStep(interpolator, lastStep);
             }
             stepStart = nextStep;
@@ -286,225 +321,33 @@
 
                 // some events handler has triggered changes that
                 // invalidate the derivatives, we need to restart from scratch
-                stopTime = start(previousF.length, stepSize, manager, stepStart, y);
-                if (Double.isNaN(previousT[0])) {
-                    return stopTime;
-                }
-                stepStart = previousT[0];
-
-                // convert to Nordsieck representation
-                scaled    = convertToNordsieckLow();
-                nordsieck = convertToNordsieckHigh(scaled);
+                start(stepStart, y, t);
                 interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
 
             }
 
+            if (! lastStep) {
+                // in some rare cases we may get here with stepSize = 0, for example
+                // when an event occurs at integration start, reducing the first step
+                // to zero; we have to reset the step to some safe non zero value
+                stepSize = filterStep(stepSize, forward, true);
+
+                // 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);
+                interpolator.rescale(hNew);
+            }
+
         }
 
-        stopTime  = stepStart;
+        final double stopTime  = stepStart;
         stepStart = Double.NaN;
         stepSize  = Double.NaN;
         return stopTime;
 
     }
 
-    /** Convert the multistep representation after a restart to Nordsieck representation.
-     * @return first scaled derivative
-     */
-    private double[] convertToNordsieckLow() {
-
-        final double[] f0 = previousF[0];
-        final double[] scaled = new double[f0.length];
-        for (int j = 0; j < f0.length; ++j) {
-            scaled[j] = stepSize * f0[j];
-        }
-        return scaled;
-
-    }
-
-    /** Convert the multistep representation after a restart to Nordsieck representation.
-     * @param scaled first scaled derivative
-     * @return Nordsieck matrix of the higher scaled derivatives
-     */
-    private RealMatrix convertToNordsieckHigh(final double[] scaled) {
-
-        final double[] f0 = previousF[0];
-        final double[][] multistep = new double[coefficients.msToN.getColumnDimension()][f0.length];
-        for (int i = 0; i < multistep.length; ++i) {
-            final double[] msI = multistep[i];
-            final double[] fI  = previousF[i + 1];
-            for (int j = 0; j < f0.length; ++j) {
-                msI[j] = stepSize * fI[j] - scaled[j];
-            }
-        }
-
-        return coefficients.msToN.multiply(new Array2DRowRealMatrix(multistep, false));
-
-    }
-
-    /** Updater for Nordsieck vector. */
-    private static class NordsieckUpdater extends DefaultRealMatrixChangingVisitor {
-
-        /** Scaled first derivative at step start. */
-        private final double[] start;
-
-        /** Scaled first derivative at step end. */
-        private final double[] end;
-
-        /** Update coefficients. */
-        private final double[] c1;
-
-        /** Simple constructor.
-         * @param start scaled first derivative at step start
-         * @param end scaled first derivative at step end
-         * @param c1 update coefficients
-         */
-        public NordsieckUpdater(final double[] start, final double[] end,
-                                final double[] c1) {
-            this.start = start;
-            this.end   = end;
-            this.c1    = c1;
-        }
-
-       /** {@inheritDoc} */
-        @Override
-        public double visit(int row, int column, double value)
-            throws MatrixVisitorException {
-            return value + c1[row] * (start[column] - end[column]);
-        }
-
-    }
-
-    /** Cache for already computed coefficients. */
-    private static class CachedCoefficients {
-
-        /** Transformer between multistep and Nordsieck representations. */
-        private final RealMatrix msToN;
-
-        /** Update coefficients of the higher order derivatives wrt y'', y''' ... */
-        private final RealMatrix msUpdate;
-
-        /** Update coefficients of the higher order derivatives wrt y'. */
-        private final double[] c1;
-
-        /** Simple constructor.
-         * @param order order of the method (must be greater than 1: due to
-         * an implementation limitation the order 1 method is not supported)
-         */
-        public CachedCoefficients(int order) {
-
-            // compute exact coefficients
-            FieldMatrix<BigFraction> bigNtoMS = buildP(order);
-            FieldMatrix<BigFraction> bigMStoN =
-                new FieldLUDecompositionImpl<BigFraction>(bigNtoMS).getSolver().getInverse();
-            BigFraction[] u = new BigFraction[order - 1];
-            Arrays.fill(u, BigFraction.ONE);
-            BigFraction[] bigC1 = bigMStoN.operate(u);
-
-            // update coefficients are computed by combining transform from
-            // Nordsieck to multistep, then shifting rows to represent step advance
-            // then applying inverse transform
-            BigFraction[][] shiftedP = bigNtoMS.getData();
-            for (int i = shiftedP.length - 1; i > 0; --i) {
-                // shift rows
-                shiftedP[i] = shiftedP[i - 1];
-            }
-            shiftedP[0] = new BigFraction[order - 1];
-            Arrays.fill(shiftedP[0], BigFraction.ZERO);
-            FieldMatrix<BigFraction> bigMSupdate =
-                bigMStoN.multiply(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
-
-            // convert coefficients to double
-            msToN    = MatrixUtils.bigFractionMatrixToRealMatrix(bigMStoN);
-            msUpdate = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
-            c1       = new double[order - 1];
-            for (int i = 0; i < order - 1; ++i) {
-                c1[i] = bigC1[i].doubleValue();
-            }
-
-        }
-
-        /** Build the P matrix transforming multistep to Nordsieck.
-         * <p>
-         * <p>
-         * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
-         * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
-         * The two representations share their two first components y(k) and
-         * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
-         * <pre>
-         * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
-         * </pre>
-         * </p>
-         * @param order order of the method (must be strictly positive)
-         * @return P matrix
-         */
-        private static FieldMatrix<BigFraction> buildP(final int order) {
-
-            final BigFraction[][] pData = new BigFraction[order - 1][order - 1];
-
-            for (int i = 0; i < pData.length; ++i) {
-                // build the P matrix elements from Taylor series formulas
-                final BigFraction[] pI = pData[i];
-                final int factor = -(i + 1);
-                int aj = factor;
-                for (int j = 0; j < pI.length; ++j) {
-                    pI[j] = new BigFraction(aj * (j + 2));
-                    aj *= factor;
-                }
-            }
-
-            return new Array2DRowFieldMatrix<BigFraction>(pData, false);
-
-        }
-
-    }
-
-    /** Serialize the instance.
-     * @param oos stream where object should be written
-     * @throws IOException if object cannot be written to stream
-     */
-    private void writeObject(ObjectOutputStream oos)
-        throws IOException {
-        oos.defaultWriteObject();
-        oos.writeInt(coefficients.msToN.getRowDimension() + 1);
-    }
-
-    /** Deserialize the instance.
-     * @param ois stream from which the object should be read
-     * @throws ClassNotFoundException if a class in the stream cannot be found
-     * @throws IOException if object cannot be read from the stream
-     */
-    private void readObject(ObjectInputStream ois)
-      throws ClassNotFoundException, IOException {
-        try {
-
-            ois.defaultReadObject();
-            final int order = ois.readInt();
-
-            final Class<AdamsBashforthIntegrator> cl = AdamsBashforthIntegrator.class;
-            final Field f = cl.getDeclaredField("coefficients");
-            f.setAccessible(true);
-
-            // cache the coefficients for each order, to avoid recomputing them
-            synchronized(cache) {
-                CachedCoefficients coeff = cache.get(order);
-                if (coeff == null) {
-                    coeff = new CachedCoefficients(order);
-                    cache.put(order, coeff);
-                }
-                f.set(this, coeff);
-            }
-
-        } catch (NoSuchFieldException nsfe) {
-            IOException ioe = new IOException();
-            ioe.initCause(nsfe);
-            throw ioe;
-        } catch (IllegalAccessException iae) {
-            IOException ioe = new IOException();
-            ioe.initCause(iae);
-            throw ioe;
-        }
-
-    }
-
 }

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Sat Jun 20 18:15:54 2009
@@ -17,25 +17,11 @@
 
 package org.apache.commons.math.ode.nonstiff;
 
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import java.io.ObjectOutputStream;
 import java.io.Serializable;
-import java.lang.reflect.Field;
 import java.util.Arrays;
-import java.util.HashMap;
-import java.util.Map;
 
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.fraction.BigFraction;
-import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
-import org.apache.commons.math.linear.FieldLUDecompositionImpl;
-import org.apache.commons.math.linear.FieldMatrix;
-import org.apache.commons.math.linear.Array2DRowFieldMatrix;
-import org.apache.commons.math.linear.MatrixUtils;
 import org.apache.commons.math.linear.MatrixVisitorException;
 import org.apache.commons.math.linear.RealMatrix;
-import org.apache.commons.math.linear.Array2DRowRealMatrix;
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
 import org.apache.commons.math.ode.DerivativeException;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
@@ -177,49 +163,51 @@
      */
     private static final long serialVersionUID = 3624292432281962886L;
 
-    /** Cache for already computed coefficients. */
-    private static final Map<Integer, CachedCoefficients> cache =
-        new HashMap<Integer, CachedCoefficients>();
-
-    /** Coefficients of the method. */
-    private final transient CachedCoefficients coefficients;
-
-    /** 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 greater than 1: due to
      * an implementation limitation the order 1 method is not supported)
-     * @param step integration step size
+     * @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 AdamsMoultonIntegrator(final int order, final double step)
+    public AdamsMoultonIntegrator(final int order,
+                                  final double minStep, final double maxStep,
+                                  final double scalAbsoluteTolerance,
+                                  final double scalRelativeTolerance)
         throws IllegalArgumentException {
+        super("Adams-Moulton", order, order, minStep, maxStep,
+              scalAbsoluteTolerance, scalRelativeTolerance);
+    }
 
-        super("Adams-Moulton", order);
-        if (order <= 1) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  "{0} is supported only for orders 2 or more",
-                  getName());
-        }
-
-        // cache the coefficients for each order, to avoid recomputing them
-        synchronized(cache) {
-            CachedCoefficients coeff = cache.get(order);
-            if (coeff == null) {
-                coeff = new CachedCoefficients(order);
-                cache.put(order, coeff);
-            }
-            coefficients = coeff;
-        }
-
-        this.step = Math.abs(step);
-
+    /**
+     * Build an Adams-Moulton integrator with the given order and step size.
+     * @param order order of the method (must be greater than 1: due to
+     * an implementation limitation the order 1 method is not supported)
+     * @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 AdamsMoultonIntegrator(final int order,
+                                  final double minStep, final double maxStep,
+                                  final double[] vecAbsoluteTolerance,
+                                  final double[] vecRelativeTolerance)
+        throws IllegalArgumentException {
+        super("Adams-Moulton", order, order, minStep, maxStep,
+              vecAbsoluteTolerance, vecRelativeTolerance);
     }
       
     
     /** {@inheritDoc} */
+    @Override
     public double integrate(final FirstOrderDifferentialEquations equations,
                             final double t0, final double[] y0,
                             final double t, final double[] y)
@@ -235,6 +223,7 @@
         if (y != y0) {
             System.arraycopy(y0, 0, y, 0, n);
         }
+        final double[] yDot = new double[y0.length];
         final double[] yTmp = new double[y0.length];
 
         // set up two interpolators sharing the integrator arrays
@@ -244,34 +233,27 @@
         interpolatorTmp.reinitialize(yTmp, forward);
 
         // set up integration control objects
-        stepStart = t0;
-        stepSize  = forward ? step : -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, stepStart, y);
-        if (Double.isNaN(previousT[0])) {
-            return stopTime;
-        }
-        stepStart = previousT[0];
-        double hNew = 0;
 
-        // convert to Nordsieck representation
-        double[]   scaled    = convertToNordsieckLow();
-        RealMatrix nordsieck = convertToNordsieckHigh(scaled);
+        // compute the initial Nordsieck vector using the configured starter integrator
+        start(t0, y, t);
         interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
         interpolator.storeTime(stepStart);
 
+        double hNew = stepSize;
+        interpolator.rescale(hNew);
+        
         boolean lastStep = false;
         while (!lastStep) {
 
             // shift all data
             interpolator.shift();
 
-            hNew  = forward ? step : -step;
+            double error = 0;
             for (boolean loop = true; loop;) {
 
                 stepSize = hNew;
@@ -282,51 +264,60 @@
                 System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
 
                 // evaluate a first estimate of the derivative (first E in the PECE sequence)
-                final double[] f0 = previousF[0];
-                previousT[0] = stepEnd;
-                computeDerivatives(stepEnd, yTmp, f0);
+                computeDerivatives(stepEnd, yTmp, yDot);
 
                 // update Nordsieck vector
-                final RealMatrix nordsieckTmp = coefficients.msUpdate.multiply(nordsieck);
                 final double[] predictedScaled = new double[y0.length];
                 for (int j = 0; j < y0.length; ++j) {
-                    predictedScaled[j] = stepSize * f0[j];
+                    predictedScaled[j] = stepSize * yDot[j];
                 }
-                nordsieckTmp.walkInOptimizedOrder(new NordsieckUpdater(scaled, predictedScaled, coefficients.c1));
+                final RealMatrix nordsieckTmp =
+                    transformer.updateHighOrderDerivativesPhase1(nordsieck);
+                transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
 
                 // apply correction (C in the PECE sequence)
-                nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
+                error = nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
 
-                // evaluate a final estimate of the derivative (second E in the PECE sequence)
-                computeDerivatives(stepEnd, yTmp, f0);
+                if (error <= 1.0) {
 
-                // update Nordsieck vector
-                final double[] correctedScaled = new double[y0.length];
-                for (int j = 0; j < y0.length; ++j) {
-                    correctedScaled[j] = stepSize * f0[j];
-                }
-                nordsieckTmp.walkInOptimizedOrder(new NordsieckUpdater(predictedScaled, correctedScaled, coefficients.c1));
+                    // evaluate a final estimate of the derivative (second E in the PECE sequence)
+                    computeDerivatives(stepEnd, yTmp, yDot);
 
-                // discrete events handling
-                interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
-                interpolatorTmp.storeTime(stepStart);
-                interpolatorTmp.shift();
-                interpolatorTmp.storeTime(stepEnd);
-                if (manager.evaluateStep(interpolatorTmp)) {
-                    final double dt = manager.getEventTime() - stepStart;
-                    if (Math.abs(dt) <= Math.ulp(stepStart)) {
-                        // rejecting the step would lead to a too small next step, we accept it
-                        loop = false;
+                    // update Nordsieck vector
+                    final double[] correctedScaled = new double[y0.length];
+                    for (int j = 0; j < y0.length; ++j) {
+                        correctedScaled[j] = stepSize * yDot[j];
+                    }
+                    transformer.updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
+
+                    // discrete events handling
+                    interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
+                    interpolatorTmp.storeTime(stepStart);
+                    interpolatorTmp.shift();
+                    interpolatorTmp.storeTime(stepEnd);
+                    if (manager.evaluateStep(interpolatorTmp)) {
+                        final double dt = manager.getEventTime() - stepStart;
+                        if (Math.abs(dt) <= Math.ulp(stepStart)) {
+                            // rejecting the step would lead to a too small next step, we accept it
+                            loop = false;
+                        } else {
+                            // reject the step to match exactly the next switch time
+                            hNew = dt;
+                            interpolator.rescale(hNew);
+                        }
                     } else {
-                        // reject the step to match exactly the next switch time
-                        hNew = dt;
+                        // accept the step
+                        scaled    = correctedScaled;
+                        nordsieck = nordsieckTmp;
+                        interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+                        loop = false;
                     }
+
                 } else {
-                    // accept the step
-                    scaled    = correctedScaled;
-                    nordsieck = nordsieckTmp;
-                    interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
-                    loop = false;
+                    // 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);
                 }
 
             }
@@ -339,8 +330,8 @@
             lastStep = manager.stop();
 
             // provide the step data to the step handler
-            interpolator.setInterpolatedTime(nextStep);
             for (StepHandler handler : stepHandlers) {
+                interpolator.setInterpolatedTime(nextStep);
                 handler.handleStep(interpolator, lastStep);
             }
             stepStart = nextStep;
@@ -349,63 +340,35 @@
 
                 // some events handler has triggered changes that
                 // invalidate the derivatives, we need to restart from scratch
-                stopTime = start(previousF.length, stepSize, manager, stepStart, y);
-                if (Double.isNaN(previousT[0])) {
-                    return stopTime;
-                }
-                stepStart = previousT[0];
-
-                // convert to Nordsieck representation
-                scaled    = convertToNordsieckLow();
-                nordsieck = convertToNordsieckHigh(scaled);
+                start(stepStart, y, t);
                 interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
-                interpolator.storeTime(stepStart);
 
             }
 
+            if (! lastStep) {
+                // in some rare cases we may get here with stepSize = 0, for example
+                // when an event occurs at integration start, reducing the first step
+                // to zero; we have to reset the step to some safe non zero value
+                stepSize = filterStep(stepSize, forward, true);
+
+                // 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);
+                interpolator.rescale(hNew);
+            }
+
         }
 
-        stopTime  = stepStart;
+        final double stopTime  = stepStart;
         stepStart = Double.NaN;
         stepSize  = Double.NaN;
         return stopTime;
 
     }
 
-    /** Convert the multistep representation after a restart to Nordsieck representation.
-     * @return first scaled derivative
-     */
-    private double[] convertToNordsieckLow() {
-
-        final double[] f0 = previousF[0];
-        final double[] scaled = new double[f0.length];
-        for (int j = 0; j < f0.length; ++j) {
-            scaled[j] = stepSize * f0[j];
-        }
-        return scaled;
-
-    }
-
-    /** Convert the multistep representation after a restart to Nordsieck representation.
-     * @param scaled first scaled derivative
-     * @return Nordsieck matrix of the higher scaled derivatives
-     */
-    private RealMatrix convertToNordsieckHigh(final double[] scaled) {
-
-        final double[] f0 = previousF[0];
-        final double[][] multistep = new double[coefficients.msToN.getColumnDimension()][f0.length];
-        for (int i = 0; i < multistep.length; ++i) {
-            final double[] msI = multistep[i];
-            final double[] fI  = previousF[i + 1];
-            for (int j = 0; j < f0.length; ++j) {
-                msI[j] = stepSize * fI[j] - scaled[j];
-            }
-        }
-
-        return coefficients.msToN.multiply(new Array2DRowRealMatrix(multistep, false));
-
-    }
-
     /** Corrector for current state in Adams-Moulton method.
      * <p>
      * This visitor implements the Taylor series formula:
@@ -414,7 +377,7 @@
      * </pre>
      * </p>
      */
-    private static class Corrector implements RealMatrixPreservingVisitor {
+    private class Corrector implements RealMatrixPreservingVisitor {
 
         /** Previous state. */
         private final double[] previous;
@@ -422,214 +385,65 @@
         /** Current scaled first derivative. */
         private final double[] scaled;
 
-        /** Placeholder where to put the recomputed current state. */
-        private final double[] corrected;
+        /** Current state before correction. */
+        private final double[] before;
+
+        /** Current state after correction. */
+        private final double[] after;
 
         /** Simple constructor.
          * @param previous previous state
          * @param scaled current scaled first derivative
-         * @param corrected placeholder where to put the corrected current state
+         * @param state state to correct (will be overwritten after visit)
          */
-        public Corrector(final double[] previous, final double[] scaled, final double[] corrected) {
+        public Corrector(final double[] previous, final double[] scaled, final double[] state) {
             this.previous = previous;
-            this.scaled    = scaled;
-            this.corrected  = corrected;
+            this.scaled   = scaled;
+            this.after    = state;
+            this.before   = state.clone();
         }
 
         /** {@inheritDoc} */
         public void start(int rows, int columns,
                           int startRow, int endRow, int startColumn, int endColumn) {
-            Arrays.fill(corrected, 0.0);
+            Arrays.fill(after, 0.0);
         }
 
         /** {@inheritDoc} */
         public void visit(int row, int column, double value)
             throws MatrixVisitorException {
             if ((row & 0x1) == 0) {
-                corrected[column] -= value;
+                after[column] -= value;
             } else {
-                corrected[column] += value;
-            }
-        }
-
-        /** {@inheritDoc} */
-        public double end() {
-            for (int i = 0; i < corrected.length; ++i) {
-                corrected[i] += previous[i] + scaled[i];
+                after[column] += value;
             }
-            return 0;
-        }
-    }
-
-    /** Updater for Nordsieck vector. */
-    private static class NordsieckUpdater extends DefaultRealMatrixChangingVisitor {
-
-        /** Scaled first derivative at step start. */
-        private final double[] start;
-
-        /** Scaled first derivative at step end. */
-        private final double[] end;
-
-        /** Update coefficients. */
-        private final double[] c1;
-
-        /** Simple constructor.
-         * @param start scaled first derivative at step start
-         * @param end scaled first derivative at step end
-         * @param c1 update coefficients
-         */
-        public NordsieckUpdater(final double[] start, final double[] end,
-                                final double[] c1) {
-            this.start = start;
-            this.end   = end;
-            this.c1    = c1;
-        }
-
-       /** {@inheritDoc} */
-        @Override
-        public double visit(int row, int column, double value)
-            throws MatrixVisitorException {
-            return value + c1[row] * (start[column] - end[column]);
         }
 
-    }
-
-    /** Cache for already computed coefficients. 
-     * @param <impements>*/
-    private static class CachedCoefficients implements Serializable {
-
         /**
-         * Serialization UID
+         * End visiting te Nordsieck vector.
+         * <p>The correction is used to control stepsize. So its amplitude is
+         * considered to be an error, which must be normalized according to
+         * error control settings. If the normalized value is greater than 1,
+         * the correction was too large and the step must be rejected.</p>
+         * @return the normalized correction, if greater than 1, the step
+         * must be rejected
          */
-        private static final long serialVersionUID = -8464316300182136812L;
-
-        /** Transformer between multistep and Nordsieck representations. */
-        private final RealMatrix msToN;
-
-        /** Update coefficients of the higher order derivatives wrt y'', y''' ... */
-        private final RealMatrix msUpdate;
-
-        /** Update coefficients of the higher order derivatives wrt y'. */
-        private final double[] c1;
-
-        /** Simple constructor.
-         * @param order order of the method (must be greater than 1: due to
-         * an implementation limitation the order 1 method is not supported)
-         */
-        public CachedCoefficients(int order) {
-
-            // compute exact coefficients
-            FieldMatrix<BigFraction> bigNtoMS = buildP(order);
-            FieldMatrix<BigFraction> bigMStoN =
-                new FieldLUDecompositionImpl<BigFraction>(bigNtoMS).getSolver().getInverse();
-            BigFraction[] u = new BigFraction[order - 1];
-            Arrays.fill(u, BigFraction.ONE);
-            BigFraction[] bigC1 = bigMStoN.operate(u);
-
-            // update coefficients are computed by combining transform from
-            // Nordsieck to multistep, then shifting rows to represent step advance
-            // then applying inverse transform
-            BigFraction[][] shiftedP = bigNtoMS.getData();
-            for (int i = shiftedP.length - 1; i > 0; --i) {
-                // shift rows
-                shiftedP[i] = shiftedP[i - 1];
-            }
-            shiftedP[0] = new BigFraction[order - 1];
-            Arrays.fill(shiftedP[0], BigFraction.ZERO);
-            FieldMatrix<BigFraction> bigMSupdate =
-                bigMStoN.multiply(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
-
-            // convert coefficients to double
-            msToN    = MatrixUtils.bigFractionMatrixToRealMatrix(bigMStoN);
-            msUpdate = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
-            c1       = new double[order - 1];
-            for (int i = 0; i < order - 1; ++i) {
-                c1[i] = bigC1[i].doubleValue();
-            }
-
-        }
-
-        /** Build the P matrix transforming multistep to Nordsieck.
-         * <p>
-         * <p>
-         * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
-         * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
-         * The two representations share their two first components y(k) and
-         * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
-         * <pre>
-         * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
-         * </pre>
-         * </p>
-         * @param order order of the method (must be strictly positive)
-         * @return P matrix
-         */
-        private static FieldMatrix<BigFraction> buildP(final int order) {
-
-            final BigFraction[][] pData = new BigFraction[order - 1][order - 1];
+        public double end() {
 
-            for (int i = 0; i < pData.length; ++i) {
-                // build the P matrix elements from Taylor series formulas
-                final BigFraction[] pI = pData[i];
-                final int factor = -(i + 1);
-                int aj = factor;
-                for (int j = 0; j < pI.length; ++j) {
-                    pI[j] = new BigFraction(aj * (j + 2));
-                    aj *= factor;
-                }
+            double error = 0;
+            for (int i = 0; i < after.length; ++i) {
+                after[i] += previous[i] + scaled[i];
+                final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i]));
+                final double tol = (vecAbsoluteTolerance == null) ?
+                                   (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
+                                   (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
+                final double ratio  = (after[i] - before[i]) / tol;
+                error += ratio * ratio;
             }
 
-            return new Array2DRowFieldMatrix<BigFraction>(pData, false);
+            return Math.sqrt(error / after.length);
 
         }
-
-    }
-
-    /** Serialize the instance.
-     * @param oos stream where object should be written
-     * @throws IOException if object cannot be written to stream
-     */
-    private void writeObject(ObjectOutputStream oos)
-        throws IOException {
-        oos.defaultWriteObject();
-        oos.writeInt(coefficients.msToN.getRowDimension() + 1);
-    }
-
-    /** Deserialize the instance.
-     * @param ois stream from which the object should be read
-     * @throws ClassNotFoundException if a class in the stream cannot be found
-     * @throws IOException if object cannot be read from the stream
-     */
-    private void readObject(ObjectInputStream ois)
-      throws ClassNotFoundException, IOException {
-        try {
-
-            ois.defaultReadObject();
-            final int order = ois.readInt();
-
-            final Class<AdamsMoultonIntegrator> cl = AdamsMoultonIntegrator.class;
-            final Field f = cl.getDeclaredField("coefficients");
-            f.setAccessible(true);
-
-            // cache the coefficients for each order, to avoid recomputing them
-            synchronized(cache) {
-                CachedCoefficients coeff = cache.get(order);
-                if (coeff == null) {
-                    coeff = new CachedCoefficients(order);
-                    cache.put(order, coeff);
-                }
-                f.set(this, coeff);
-            }
-
-        } catch (NoSuchFieldException nsfe) {
-            IOException ioe = new IOException();
-            ioe.initCause(nsfe);
-            throw ioe;
-        } catch (IllegalAccessException iae) {
-            IOException ioe = new IOException();
-            ioe.initCause(iae);
-            throw ioe;
-        }
-
     }
 
 }

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java Sat Jun 20 18:15:54 2009
@@ -22,8 +22,9 @@
 import java.io.ObjectOutput;
 import java.util.Arrays;
 
-import org.apache.commons.math.linear.RealMatrix;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
+import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
+import org.apache.commons.math.linear.RealMatrix;
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
 
 /**
@@ -93,7 +94,7 @@
         return new NordsieckStepInterpolator(this);
     }
 
-    /** Reinitialize the instance
+    /** Reinitialize the instance.
      * <p>Beware that all arrays <em>must</em> be references to integrator
      * arrays, in order to ensure proper update without copy.</p>
      * @param y reference to the integrator array holding the state at
@@ -105,7 +106,7 @@
         super.reinitialize(y, forward);
     }
 
-    /** Reinitialize the instance
+    /** Reinitialize the instance.
      * <p>Beware that all arrays <em>must</em> be references to integrator
      * arrays, in order to ensure proper update without copy.</p>
      * @param referenceTime time at which all arrays are defined
@@ -127,6 +128,20 @@
 
     }
 
+    /** Rescale the instance.
+     * <p>Since the scaled and Nordiseck arrays are shared with the caller,
+     * this method has the side effect of rescaling this arrays in the caller too.</p>
+     * @param scalingH new step size to use in the scaled and nordsieck arrays
+     */
+    public void rescale(final double scalingH) {
+        final double ratio = scalingH / this.scalingH;
+        for (int i = 0; i < scaled.length; ++i) {
+            scaled[i] *= ratio;
+        }
+        nordsieck.walkInOptimizedOrder(new Rescaler(ratio));
+        this.scalingH = scalingH;
+    }
+
     /** {@inheritDoc} */
     @Override
     protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
@@ -151,7 +166,7 @@
          * @param theta normalized interpolation abscissa within the step
          */
         public StateEstimator(final double scale, final double theta) {
-            this.scale  = scale;
+            this.scale = scale;
             lowPower   = theta;
             highPowers = new double[nordsieck.getRowDimension()];
             double thetaN = theta;
@@ -187,6 +202,32 @@
 
     }
 
+    /** Visitor rescaling the Nordsieck vector. */
+    private class Rescaler extends DefaultRealMatrixChangingVisitor {
+
+        /** Powers of the rescaling ratio. */
+        private final double[] powers;
+
+        /** Simple constructor.
+         * @param ratio rescaling ratio
+         */
+        public Rescaler(final double ratio) {
+            powers = new double[nordsieck.getRowDimension()];
+            double f = ratio;
+            for (int i = 0; i < powers.length; ++i) {
+                f *= ratio;
+                powers[i] = f;
+            }
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public double visit(final int row, final int column, final double value) {
+            return value * powers[row];
+        }
+        
+    }
+
     /** {@inheritDoc} */
     @Override
     public void writeExternal(final ObjectOutput out)

Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Sat Jun 20 18:15:54 2009
@@ -29,7 +29,6 @@
 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;
 import org.junit.Test;
 
 public class AdamsBashforthIntegratorTest {
@@ -37,70 +36,67 @@
     @Test(expected=IntegratorException.class)
     public void dimensionCheck() throws DerivativeException, IntegratorException {
         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]);
+        FirstOrderIntegrator integ =
+            new AdamsBashforthIntegrator(3, 0.0, 1.0, 1.0e-10, 1.0e-10);
+        integ.integrate(pb,
+                        0.0, new double[pb.getDimension()+10],
+                        1.0, new double[pb.getDimension()+10]);
     }
 
-    @Test
-    public void decreasingSteps() 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-3 * 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) && !(pb instanceof TestProblem4) && !(pb instanceof TestProblem6)) {
-                    assertTrue(error <= Math.abs(1.05 * previousError));
-                }
-                previousError = error;
+    @Test(expected=IntegratorException.class)
+    public void testMinStep() throws DerivativeException, IntegratorException {
 
-            }
-
-        }
+          TestProblem1 pb = new TestProblem1();
+          double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
+          double maxStep = pb.getFinalTime() - pb.getInitialTime();
+          double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
+          double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
+
+          FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, minStep, maxStep,
+                                                                    vecAbsoluteTolerance,
+                                                                    vecRelativeTolerance);
+          TestProblemHandler handler = new TestProblemHandler(pb, integ);
+          integ.addStepHandler(handler);
+          integ.integrate(pb,
+                          pb.getInitialTime(), pb.getInitialState(),
+                          pb.getFinalTime(), new double[pb.getDimension()]);
 
     }
 
     @Test
-    public void smallStep() throws DerivativeException, IntegratorException {
-
-        TestProblem1 pb  = new TestProblem1();
-        double range = pb.getFinalTime() - pb.getInitialTime();
-        double step = range * 0.001;
+    public void testIncreasingTolerance()
+        throws DerivativeException, IntegratorException {
+
+        int previousCalls = Integer.MAX_VALUE;
+        for (int i = -12; i < -2; ++i) {
+            TestProblem1 pb = new TestProblem1();
+            double minStep = 0;
+            double maxStep = pb.getFinalTime() - pb.getInitialTime();
+            double scalAbsoluteTolerance = Math.pow(10.0, i);
+            double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+            FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, minStep, maxStep,
+                                                                      scalAbsoluteTolerance,
+                                                                      scalRelativeTolerance);
+            TestProblemHandler handler = new TestProblemHandler(pb, integ);
+            integ.addStepHandler(handler);
+            integ.integrate(pb,
+                            pb.getInitialTime(), pb.getInitialState(),
+                            pb.getFinalTime(), new double[pb.getDimension()]);
 
-        AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, step);
-        integ.setStarterIntegrator(new DormandPrince853Integrator(0, range, 1.0e-12, 1.0e-12));
-        TestProblemHandler handler = new TestProblemHandler(pb, integ);
-        integ.addStepHandler(handler);
-        integ.integrate(pb,
-                        pb.getInitialTime(), pb.getInitialState(),
-                        pb.getFinalTime(), new double[pb.getDimension()]);
+            // the 28 and 42 factors are only valid for this test
+            // and has been obtained from trial and error
+            // there is no general relation between local and global errors
+            assertTrue(handler.getMaximalValueError() > (28.0 * scalAbsoluteTolerance));
+            assertTrue(handler.getMaximalValueError() < (42.0 * scalAbsoluteTolerance));
+            assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
+
+            int calls = pb.getCalls();
+            assertEquals(integ.getEvaluations(), calls);
+            assertTrue(calls <= previousCalls);
+            previousCalls = calls;
 
-        assertTrue(handler.getLastError() < 2.0e-9);
-        assertTrue(handler.getMaximalValueError() < 9.0e-9);
-        assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
-        assertEquals("Adams-Bashforth", integ.getName());
-        assertTrue(integ.getEvaluations() > 1000);
-        assertEquals(Integer.MAX_VALUE, integ.getMaxEvaluations());
+        }
 
     }
 
@@ -109,13 +105,11 @@
 
         TestProblem1 pb  = new TestProblem1();
         double range = pb.getFinalTime() - pb.getInitialTime();
-        double step = range * 0.001;
 
-        AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, step);
-        integ.setStarterIntegrator(new DormandPrince853Integrator(0, range, 1.0e-12, 1.0e-12));
+        AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, 0, range, 1.0e-12, 1.0e-12);
         TestProblemHandler handler = new TestProblemHandler(pb, integ);
         integ.addStepHandler(handler);
-        integ.setMaxEvaluations(1000);
+        integ.setMaxEvaluations(650);
         integ.integrate(pb,
                         pb.getInitialTime(), pb.getInitialState(),
                         pb.getFinalTime(), new double[pb.getDimension()]);
@@ -123,38 +117,19 @@
     }
 
     @Test
-    public void bigStep() 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.06);
-        assertTrue(handler.getMaximalValueError() > 0.06);
-        assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
-
-    }
-
-    @Test
     public void backward() throws DerivativeException, IntegratorException {
 
         TestProblem5 pb = new TestProblem5();
-        double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+        double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
 
-        FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, step);
+        FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, 0, range, 1.0e-12, 1.0e-12);
         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() < 8.0e-11);
-        assertTrue(handler.getMaximalValueError() < 8.0e-11);
+        assertTrue(handler.getLastError() < 1.0e-8);
+        assertTrue(handler.getMaximalValueError() < 1.0e-8);
         assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
         assertEquals("Adams-Bashforth", integ.getName());
     }
@@ -162,20 +137,19 @@
     @Test
     public void polynomial() throws DerivativeException, IntegratorException {
         TestProblem6 pb = new TestProblem6();
-        double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.02;
+        double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
 
         for (int order = 2; order < 9; ++order) {
-            AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(order, step);
-            integ.setStarterIntegrator(new DormandPrince853Integrator(1.0e-3 * step, 1.0e3 * step,
-                                                                      1.0e-5, 1.0e-5));
+            AdamsBashforthIntegrator integ =
+                new AdamsBashforthIntegrator(order, 1.0e-6 * range, 0.1 * range, 1.0e-10, 1.0e-10);
             TestProblemHandler handler = new TestProblemHandler(pb, integ);
             integ.addStepHandler(handler);
             integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
                             pb.getFinalTime(), new double[pb.getDimension()]);
             if (order < 5) {
-                assertTrue(handler.getMaximalValueError() > 1.0e-5);
+                assertTrue(integ.getEvaluations() > 160);
             } else {
-                assertTrue(handler.getMaximalValueError() < 7.0e-12);
+                assertTrue(integ.getEvaluations() < 70);
             }
         }
 
@@ -187,12 +161,14 @@
                IOException, ClassNotFoundException {
 
         TestProblem6 pb = new TestProblem6();
-        double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.01;
+        double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
 
         ByteArrayOutputStream bos = new ByteArrayOutputStream();
         ObjectOutputStream    oos = new ObjectOutputStream(bos);
-        oos.writeObject(new AdamsBashforthIntegrator(8, step));
-        
+        oos.writeObject(new AdamsBashforthIntegrator(5, 0, range, 1.0e-12, 1.0e-12));
+        assertTrue(bos.size() > 2800);
+        assertTrue(bos.size() < 3000);
+
         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
         ObjectInputStream     ois = new ObjectInputStream(bis);
         FirstOrderIntegrator integ  = (AdamsBashforthIntegrator) ois.readObject();
@@ -201,7 +177,7 @@
         integ.addStepHandler(handler);
         integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
                         pb.getFinalTime(), new double[pb.getDimension()]);
-        assertTrue(handler.getMaximalValueError() < 7.0e-13);
+        assertTrue(handler.getMaximalValueError() < 2.0e-11);
 
     }
 



Mime
View raw message