commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject [3/4] [math] Base classes for field-based multistep integrators.
Date Mon, 28 Dec 2015 21:14:43 GMT
Base classes for field-based multistep integrators.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/ba6bd8cb
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/ba6bd8cb
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/ba6bd8cb

Branch: refs/heads/MATH_3_X
Commit: ba6bd8cbe193b86d9fa44b1c68c417ab88953e29
Parents: 4a501b6
Author: Luc Maisonobe <luc@apache.org>
Authored: Mon Dec 28 22:13:48 2015 +0100
Committer: Luc Maisonobe <luc@apache.org>
Committed: Mon Dec 28 22:13:48 2015 +0100

----------------------------------------------------------------------
 .../math3/ode/MultistepFieldIntegrator.java     | 427 +++++++++++++++++++
 .../ode/nonstiff/AdamsFieldIntegrator.java      | 146 +++++++
 .../nonstiff/AdamsFieldStepInterpolator.java    | 193 +++++++++
 .../AdamsNordsieckFieldTransformer.java         | 362 ++++++++++++++++
 .../AbstractAdamsFieldIntegratorTest.java       | 261 ++++++++++++
 5 files changed, 1389 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
new file mode 100644
index 0000000..b86fc05
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
@@ -0,0 +1,427 @@
+/*
+ * 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.math3.ode;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MathIllegalStateException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
+import org.apache.commons.math3.ode.nonstiff.DormandPrince853FieldIntegrator;
+import org.apache.commons.math3.ode.sampling.FieldStepHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+import org.apache.commons.math3.util.MathUtils;
+
+/**
+ * This class is the base class for multistep integrators for Ordinary
+ * Differential Equations.
+ * <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<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
+ * </pre></p>
+ * <p>Rather than storing several previous steps separately, this implementation uses
+ * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
+ * step (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>
+ * (we omit the k index in the notation for clarity)</p>
+ * <p>
+ * Multistep integrators with Nordsieck representation are highly sensitive to
+ * large step changes because when the step is multiplied by factor a, the
+ * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
+ * and the last components are the least accurate ones. The default max growth
+ * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
+ * </p>
+ *
+ * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthFieldIntegrator
+ * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonFieldIntegrator
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
+    extends AdaptiveStepsizeFieldIntegrator<T> {
+
+    /** First scaled derivative (h y'). */
+    protected T[] scaled;
+
+    /** Nordsieck matrix of the higher scaled derivatives.
+     * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
+     */
+    protected Array2DRowFieldMatrix<T> nordsieck;
+
+    /** Starter integrator. */
+    private FirstOrderFieldIntegrator<T> starter;
+
+    /** Number of steps of the multistep method (excluding the one being computed). */
+    private final int nSteps;
+
+    /** Stepsize control exponent. */
+    private double exp;
+
+    /** Safety factor for stepsize control. */
+    private double safety;
+
+    /** Minimal reduction factor for stepsize control. */
+    private double minReduction;
+
+    /** Maximal growth factor for stepsize control. */
+    private double maxGrowth;
+
+    /**
+     * Build a multistep integrator with the given stepsize bounds.
+     * <p>The default starter integrator is set to the {@link
+     * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
+     * some defaults settings.</p>
+     * <p>
+     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
+     * </p>
+     * @param field field to which the time and state vector elements belong
+     * @param name name of the method
+     * @param nSteps number of steps of the multistep method
+     * (excluding the one being computed)
+     * @param order order of the method
+     * @param minStep minimal step (must be positive even for backward
+     * integration), the last step can be smaller than this
+     * @param maxStep maximal step (must be positive even for backward
+     * integration)
+     * @param scalAbsoluteTolerance allowed absolute error
+     * @param scalRelativeTolerance allowed relative error
+     * @exception NumberIsTooSmallException if number of steps is smaller than 2
+     */
+    protected MultistepFieldIntegrator(final Field<T> field, final String name,
+                                       final int nSteps, final int order,
+                                       final double minStep, final double maxStep,
+                                       final double scalAbsoluteTolerance,
+                                       final double scalRelativeTolerance)
+        throws NumberIsTooSmallException {
+
+        super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+
+        if (nSteps < 2) {
+            throw new NumberIsTooSmallException(
+                  LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
+                  nSteps, 2, true);
+        }
+
+        starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
+                                                         scalAbsoluteTolerance,
+                                                         scalRelativeTolerance);
+        this.nSteps = nSteps;
+
+        exp = -1.0 / order;
+
+        // set the default values of the algorithm control parameters
+        setSafety(0.9);
+        setMinReduction(0.2);
+        setMaxGrowth(FastMath.pow(2.0, -exp));
+
+    }
+
+    /**
+     * Build a multistep integrator with the given stepsize bounds.
+     * <p>The default starter integrator is set to the {@link
+     * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
+     * some defaults settings.</p>
+     * <p>
+     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
+     * </p>
+     * @param field field to which the time and state vector elements belong
+     * @param name name of the method
+     * @param nSteps number of steps of the multistep method
+     * (excluding the one being computed)
+     * @param order order of the method
+     * @param minStep minimal step (must be positive even for backward
+     * integration), the last step can be smaller than this
+     * @param maxStep maximal step (must be positive even for backward
+     * integration)
+     * @param vecAbsoluteTolerance allowed absolute error
+     * @param vecRelativeTolerance allowed relative error
+     */
+    protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
+                                       final int order,
+                                       final double minStep, final double maxStep,
+                                       final double[] vecAbsoluteTolerance,
+                                       final double[] vecRelativeTolerance) {
+        super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+        starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
+                                                         vecAbsoluteTolerance,
+                                                         vecRelativeTolerance);
+        this.nSteps = nSteps;
+
+        exp = -1.0 / order;
+
+        // set the default values of the algorithm control parameters
+        setSafety(0.9);
+        setMinReduction(0.2);
+        setMaxGrowth(FastMath.pow(2.0, -exp));
+
+    }
+
+    /**
+     * Get the starter integrator.
+     * @return starter integrator
+     */
+    public FirstOrderFieldIntegrator<T> 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 starterIntegrator starter integrator
+     */
+    public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) {
+        this.starter = starterIntegrator;
+    }
+
+    /** Start the integration.
+     * <p>This method computes one step using the underlying starter integrator,
+     * and initializes the Nordsieck vector at step start. The starter integrator
+     * purpose is only to establish initial conditions, it does not really change
+     * time by itself. The top level multistep integrator remains in charge of
+     * handling time propagation and events handling as it will starts its own
+     * computation right from the beginning. In a sense, the starter integrator
+     * can be seen as a dummy one and so it will never trigger any user event nor
+     * call any user step handler.</p>
+     * @param equations complete set of differential equations to integrate
+     * @param initialState initial state (time, primary and secondary state vectors)
+     * @param t target time for the integration
+     * (can be set to a value smaller than <code>t0</code> for backward integration)
+     * @exception DimensionMismatchException if arrays dimension do not match equations settings
+     * @exception NumberIsTooSmallException if integration step is too small
+     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+     * @exception NoBracketingException if the location of an event cannot be bracketed
+     */
+    protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
+        throws DimensionMismatchException, NumberIsTooSmallException,
+               MaxCountExceededException, NoBracketingException {
+
+        // make sure NO user event nor user step handler is triggered,
+        // this is the task of the top level integrator, not the task
+        // of the starter integrator
+        starter.clearEventHandlers();
+        starter.clearStepHandlers();
+
+        // set up one specific step handler to extract initial Nordsieck vector
+        starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2));
+
+        // start integration, expecting a InitializationCompletedMarkerException
+        try {
+
+            starter.integrate(equations, initialState, t);
+
+            // we should not reach this step
+            throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
+
+        } catch (InitializationCompletedMarkerException icme) { // NOPMD
+            // this is the expected nominal interruption of the start integrator
+
+            // count the evaluations used by the starter
+            getEvaluationsCounter().increment(starter.getEvaluations());
+
+        }
+
+        // remove the specific step handler
+        starter.clearStepHandlers();
+
+    }
+
+    /** Initialize the high order scaled derivatives at step start.
+     * @param h step size to use for scaling
+     * @param t first steps times
+     * @param y first steps states
+     * @param yDot first steps derivatives
+     * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
+     * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
+     */
+    protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
+                                                                               final T[][] y,
+                                                                               final T[][] yDot);
+
+    /** Get the minimal reduction factor for stepsize control.
+     * @return minimal reduction factor
+     */
+    public double getMinReduction() {
+        return minReduction;
+    }
+
+    /** Set the minimal reduction factor for stepsize control.
+     * @param minReduction minimal reduction factor
+     */
+    public void setMinReduction(final double minReduction) {
+        this.minReduction = minReduction;
+    }
+
+    /** Get the maximal growth factor for stepsize control.
+     * @return maximal growth factor
+     */
+    public double getMaxGrowth() {
+        return maxGrowth;
+    }
+
+    /** Set the maximal growth factor for stepsize control.
+     * @param maxGrowth maximal growth factor
+     */
+    public void setMaxGrowth(final double maxGrowth) {
+        this.maxGrowth = maxGrowth;
+    }
+
+    /** Get the safety factor for stepsize control.
+     * @return safety factor
+     */
+    public double getSafety() {
+      return safety;
+    }
+
+    /** Set the safety factor for stepsize control.
+     * @param safety safety factor
+     */
+    public void setSafety(final double safety) {
+      this.safety = safety;
+    }
+
+    /** Get the number of steps of the multistep method (excluding the one being computed).
+     * @return number of steps of the multistep method (excluding the one being computed)
+     */
+    public int getNSteps() {
+      return nSteps;
+    }
+
+    /** Compute step grow/shrink factor according to normalized error.
+     * @param error normalized error of the current step
+     * @return grow/shrink factor for next step
+     */
+    protected T computeStepGrowShrinkFactor(final T error) {
+        return MathUtils.min(error.getField().getZero().add(maxGrowth),
+                             MathUtils.max(error.getField().getZero().add(minReduction),
+                                           error.pow(exp).multiply(safety)));
+    }
+
+    /** Specialized step handler storing the first step.
+     */
+    private class FieldNordsieckInitializer implements FieldStepHandler<T> {
+
+        /** Equation mapper. */
+        private final FieldEquationsMapper<T> mapper;
+
+        /** Steps counter. */
+        private int count;
+
+        /** Saved start. */
+        private FieldODEStateAndDerivative<T> savedStart;
+
+        /** First steps times. */
+        private final T[] t;
+
+        /** First steps states. */
+        private final T[][] y;
+
+        /** First steps derivatives. */
+        private final T[][] yDot;
+
+        /** Simple constructor.
+         * @param mapper equation mapper
+         * @param nbStartPoints number of start points (including the initial point)
+         */
+        FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) {
+            this.mapper = mapper;
+            this.count  = 0;
+            this.t      = MathArrays.buildArray(getField(), nbStartPoints);
+            this.y      = MathArrays.buildArray(getField(), nbStartPoints, -1);
+            this.yDot   = MathArrays.buildArray(getField(), nbStartPoints, -1);
+        }
+
+        /** {@inheritDoc} */
+        public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
+            throws MaxCountExceededException {
+
+
+            if (count == 0) {
+                // first step, we need to store also the point at the beginning of the step
+                final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
+                savedStart  = prev;
+                t[count]    = prev.getTime();
+                y[count]    = mapper.mapState(prev);
+                yDot[count] = mapper.mapDerivative(prev);
+            }
+
+            // store the point at the end of the step
+            ++count;
+            final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
+            t[count]    = curr.getTime();
+            y[count]    = mapper.mapState(curr);
+            yDot[count] = mapper.mapDerivative(curr);
+
+            if (count == t.length - 1) {
+
+                // this was the last point we needed, we can compute the derivatives
+                setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));
+
+                // first scaled derivative
+                scaled = MathArrays.buildArray(getField(), yDot[0].length);
+                for (int j = 0; j < scaled.length; ++j) {
+                    scaled[j] = yDot[0][j].multiply(getStepSize());
+                }
+
+                // higher order derivatives
+                nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);
+
+                // stop the integrator now that all needed steps have been handled
+                setStepStart(savedStart);
+                throw new InitializationCompletedMarkerException();
+
+            }
+
+        }
+
+        /** {@inheritDoc} */
+        public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
+            // nothing to do
+        }
+
+    }
+
+    /** Marker exception used ONLY to stop the starter integrator after first step. */
+    private static class InitializationCompletedMarkerException
+        extends RuntimeException {
+
+        /** Serializable version identifier. */
+        private static final long serialVersionUID = -1914085471038046418L;
+
+        /** Simple constructor. */
+        InitializationCompletedMarkerException() {
+            super((Throwable) null);
+        }
+
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
new file mode 100644
index 0000000..1657ca5
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
@@ -0,0 +1,146 @@
+/*
+ * 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.math3.ode.nonstiff;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.MultistepFieldIntegrator;
+
+
+/** Base class for {@link AdamsBashforthFieldIntegrator Adams-Bashforth} and
+ * {@link AdamsMoultonFieldIntegrator Adams-Moulton} integrators.
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+public abstract class AdamsFieldIntegrator<T extends RealFieldElement<T>> extends MultistepFieldIntegrator<T> {
+
+    /** Transformer. */
+    private final AdamsNordsieckFieldTransformer<T> transformer;
+
+    /**
+     * Build an Adams integrator with the given order and step control parameters.
+     * @param field field to which the time and state vector elements belong
+     * @param name name of the method
+     * @param nSteps number of steps of the method excluding the one being computed
+     * @param order order of the method
+     * @param minStep minimal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param maxStep maximal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param scalAbsoluteTolerance allowed absolute error
+     * @param scalRelativeTolerance allowed relative error
+     * @exception NumberIsTooSmallException if order is 1 or less
+     */
+    public AdamsFieldIntegrator(final Field<T> field, final String name,
+                                final int nSteps, final int order,
+                                final double minStep, final double maxStep,
+                                final double scalAbsoluteTolerance,
+                                final double scalRelativeTolerance)
+        throws NumberIsTooSmallException {
+        super(field, name, nSteps, order, minStep, maxStep,
+              scalAbsoluteTolerance, scalRelativeTolerance);
+        transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
+    }
+
+    /**
+     * Build an Adams integrator with the given order and step control parameters.
+     * @param field field to which the time and state vector elements belong
+     * @param name name of the method
+     * @param nSteps number of steps of the method excluding the one being computed
+     * @param order order of the method
+     * @param minStep minimal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param maxStep maximal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param vecAbsoluteTolerance allowed absolute error
+     * @param vecRelativeTolerance allowed relative error
+     * @exception IllegalArgumentException if order is 1 or less
+     */
+    public AdamsFieldIntegrator(final Field<T> field, final String name,
+                                final int nSteps, final int order,
+                                final double minStep, final double maxStep,
+                                final double[] vecAbsoluteTolerance,
+                                final double[] vecRelativeTolerance)
+        throws IllegalArgumentException {
+        super(field, name, nSteps, order, minStep, maxStep,
+              vecAbsoluteTolerance, vecRelativeTolerance);
+        transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public abstract FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
+                                                            final FieldODEState<T> initialState,
+                                                            final T finalTime)
+        throws NumberIsTooSmallException, DimensionMismatchException,
+               MaxCountExceededException, NoBracketingException;
+
+    /** {@inheritDoc} */
+    @Override
+    protected Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
+                                                                      final T[][] y,
+                                                                      final T[][] yDot) {
+        return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
+    }
+
+    /** Update the high order scaled derivatives for Adams integrators (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[], Array2DRowFieldMatrix)
+     */
+    public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
+        return transformer.updateHighOrderDerivativesPhase1(highOrder);
+    }
+
+    /** Update the high order scaled derivatives Adams integrators (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(Array2DRowFieldMatrix)
+     */
+    public void updateHighOrderDerivativesPhase2(final T[] start, final T[] end,
+                                                 final Array2DRowFieldMatrix<T> highOrder) {
+        transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
new file mode 100644
index 0000000..9f282c0
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
@@ -0,0 +1,193 @@
+/*
+ * 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.math3.ode.nonstiff;
+
+import java.util.Arrays;
+
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.ode.FieldEquationsMapper;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator;
+import org.apache.commons.math3.util.MathArrays;
+
+/**
+ * This class implements an interpolator for Adams integrators using Nordsieck representation.
+ *
+ * <p>This interpolator computes dense output around the current point.
+ * The interpolation equation is based on Taylor series formulas.
+ *
+ * @see AdamsBashforthFieldIntegrator
+ * @see AdamsMoultonFieldIntegrator
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+class AdamsFieldStepInterpolator<T extends RealFieldElement<T>> extends AbstractFieldStepInterpolator<T> {
+
+    /** Step size used in the first scaled derivative and Nordsieck vector. */
+    private T scalingH;
+
+    /** First scaled derivative. */
+    private final T[] scaled;
+
+    /** Nordsieck vector. */
+    private final Array2DRowFieldMatrix<T> nordsieck;
+
+    /** Simple constructor.
+     * @param stepSize step size used in the scaled and Nordsieck arrays
+     * @param referenceState reference state from which Taylor expansion are estimated
+     * @param scaled first scaled derivative
+     * @param nordsieck Nordsieck vector
+     * @param isForward integration direction indicator
+     * @param equationsMapper mapper for ODE equations primary and secondary components
+     */
+    AdamsFieldStepInterpolator(final T stepSize, final FieldODEStateAndDerivative<T> referenceState,
+                               final T[] scaled, final Array2DRowFieldMatrix<T> nordsieck,
+                               final boolean isForward, final FieldEquationsMapper<T> equationsMapper) {
+        this(stepSize, scaled, nordsieck, isForward,
+             referenceState, taylor(referenceState, referenceState.getTime().add(stepSize), stepSize, scaled, nordsieck),
+             equationsMapper);
+    }
+
+    /** Simple constructor.
+     * @param stepSize step size used in the scaled and Nordsieck arrays
+     * @param scaled first scaled derivative
+     * @param nordsieck Nordsieck vector
+     * @param isForward integration direction indicator
+     * @param globalPreviousState start of the global step
+     * @param globalCurrentState end of the global step
+     * @param equationsMapper mapper for ODE equations primary and secondary components
+     */
+    private AdamsFieldStepInterpolator(final T stepSize, final T[] scaled,
+                                       final Array2DRowFieldMatrix<T> nordsieck,
+                                       final boolean isForward,
+                                       final FieldODEStateAndDerivative<T> globalPreviousState,
+                                       final FieldODEStateAndDerivative<T> globalCurrentState,
+                                       final FieldEquationsMapper<T> equationsMapper) {
+        this(stepSize, scaled, nordsieck,
+             isForward, globalPreviousState, globalCurrentState,
+             globalPreviousState, globalCurrentState, equationsMapper);
+    }
+
+    /** Simple constructor.
+     * @param stepSize step size used in the scaled and Nordsieck arrays
+     * @param scaled first scaled derivative
+     * @param nordsieck Nordsieck vector
+     * @param isForward integration direction indicator
+     * @param globalPreviousState start of the global step
+     * @param globalCurrentState end of the global step
+     * @param softPreviousState start of the restricted step
+     * @param softCurrentState end of the restricted step
+     * @param equationsMapper mapper for ODE equations primary and secondary components
+     */
+    private AdamsFieldStepInterpolator(final T stepSize, final T[] scaled,
+                                       final Array2DRowFieldMatrix<T> nordsieck,
+                                       final boolean isForward,
+                                       final FieldODEStateAndDerivative<T> globalPreviousState,
+                                       final FieldODEStateAndDerivative<T> globalCurrentState,
+                                       final FieldODEStateAndDerivative<T> softPreviousState,
+                                       final FieldODEStateAndDerivative<T> softCurrentState,
+                                       final FieldEquationsMapper<T> equationsMapper) {
+        super(isForward, globalPreviousState, globalCurrentState,
+              softPreviousState, softCurrentState, equationsMapper);
+        this.scalingH  = stepSize;
+        this.scaled    = scaled.clone();
+        this.nordsieck = new Array2DRowFieldMatrix<T>(nordsieck.getData(), false);
+    }
+
+    /** Create a new instance.
+     * @param newForward integration direction indicator
+     * @param newGlobalPreviousState start of the global step
+     * @param newGlobalCurrentState end of the global step
+     * @param newSoftPreviousState start of the restricted step
+     * @param newSoftCurrentState end of the restricted step
+     * @param newMapper equations mapper for the all equations
+     * @return a new instance
+     */
+    protected AdamsFieldStepInterpolator<T> create(boolean newForward,
+                                                   FieldODEStateAndDerivative<T> newGlobalPreviousState,
+                                                   FieldODEStateAndDerivative<T> newGlobalCurrentState,
+                                                   FieldODEStateAndDerivative<T> newSoftPreviousState,
+                                                   FieldODEStateAndDerivative<T> newSoftCurrentState,
+                                                   FieldEquationsMapper<T> newMapper) {
+        return new AdamsFieldStepInterpolator<T>(scalingH, scaled, nordsieck,
+                                                 newForward,
+                                                 newGlobalPreviousState, newGlobalCurrentState,
+                                                 newSoftPreviousState, newSoftCurrentState,
+                                                 newMapper);
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> equationsMapper,
+                                                                                   final T time, final T theta,
+                                                                                   final T thetaH, final T oneMinusThetaH) {
+        return taylor(getPreviousState(), time, scalingH, scaled, nordsieck);
+    }
+
+    /** Estimate state by applying Taylor formula.
+     * @param referenceState reference state
+     * @param time time at which state must be estimated
+     * @param stepSize step size used in the scaled and Nordsieck arrays
+     * @param scaled first scaled derivative
+     * @param nordsieck Nordsieck vector
+     * @return estimated state
+     * @param <S> the type of the field elements
+     */
+    private static <S extends RealFieldElement<S>> FieldODEStateAndDerivative<S> taylor(final FieldODEStateAndDerivative<S> referenceState,
+                                                                                        final S time, final S stepSize,
+                                                                                        final S[] scaled,
+                                                                                        final Array2DRowFieldMatrix<S> nordsieck) {
+
+        final S x = time.subtract(referenceState.getTime());
+        final S normalizedAbscissa = x.divide(stepSize);
+
+        S[] stateVariation = MathArrays.buildArray(time.getField(), scaled.length);
+        Arrays.fill(stateVariation, time.getField().getZero());
+        S[] estimatedDerivatives = MathArrays.buildArray(time.getField(), scaled.length);
+        Arrays.fill(estimatedDerivatives, time.getField().getZero());
+
+        // apply Taylor formula from high order to low order,
+        // for the sake of numerical accuracy
+        final S[][] nData = nordsieck.getDataRef();
+        for (int i = nData.length - 1; i >= 0; --i) {
+            final int order = i + 2;
+            final S[] nDataI = nData[i];
+            final S power = normalizedAbscissa.pow(order);
+            for (int j = 0; j < nDataI.length; ++j) {
+                final S d = nDataI[j].multiply(power);
+                stateVariation[j]          = stateVariation[j].add(d);
+                estimatedDerivatives[j] = estimatedDerivatives[j].add(d.multiply(order));
+            }
+        }
+
+        S[] estimatedState = referenceState.getState();
+        for (int j = 0; j < stateVariation.length; ++j) {
+            stateVariation[j]    = stateVariation[j].add(scaled[j].multiply(normalizedAbscissa));
+            estimatedState[j] = estimatedState[j].add(stateVariation[j]);
+            estimatedDerivatives[j] =
+                estimatedDerivatives[j].add(scaled[j].multiply(normalizedAbscissa)).divide(x);
+        }
+
+        return new FieldODEStateAndDerivative<S>(time, estimatedState, estimatedDerivatives);
+
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
new file mode 100644
index 0000000..89f8393
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
@@ -0,0 +1,362 @@
+/*
+ * 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.math3.ode.nonstiff;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.linear.ArrayFieldVector;
+import org.apache.commons.math3.linear.FieldDecompositionSolver;
+import org.apache.commons.math3.linear.FieldLUDecomposition;
+import org.apache.commons.math3.linear.FieldMatrix;
+import org.apache.commons.math3.util.MathArrays;
+
+/** Transformer to Nordsieck vectors for Adams integrators.
+ * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
+ * {@link AdamsMoultonIntegrator Adams-Moulton} 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<sup>(k)</sup><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&gt;0</sub> (j+1) (-i)<sup>j</sup> s<sub>j+1</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+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
+ * the column number starting from 1:
+ * <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 matrix is simply
+ * the absolute value of matrix P.</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>
+ * 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>
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+public class AdamsNordsieckFieldTransformer<T extends RealFieldElement<T>> {
+
+    /** Cache for already computed coefficients. */
+    private static final Map<Integer,
+                         Map<Field<? extends RealFieldElement<?>>,
+                                   AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>> CACHE =
+        new HashMap<Integer, Map<Field<? extends RealFieldElement<?>>,
+                                 AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>>();
+
+    /** Field to which the time and state vector elements belong. */
+    private final Field<T> field;
+
+    /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
+    private final Array2DRowFieldMatrix<T> update;
+
+    /** Update coefficients of the higher order derivatives wrt y'. */
+    private final T[] c1;
+
+    /** Simple constructor.
+     * @param field field to which the time and state vector elements belong
+     * @param n number of steps of the multistep method
+     * (excluding the one being computed)
+     */
+    private AdamsNordsieckFieldTransformer(final Field<T> field, final int n) {
+
+        this.field = field;
+        final int rows = n - 1;
+
+        // compute coefficients
+        FieldMatrix<T> bigP = buildP(rows);
+        FieldDecompositionSolver<T> pSolver =
+            new FieldLUDecomposition<T>(bigP).getSolver();
+
+        T[] u = MathArrays.buildArray(field, rows);
+        Arrays.fill(u, field.getOne());
+        c1 = pSolver.solve(new ArrayFieldVector<T>(u, false)).toArray();
+
+        // update coefficients are computed by combining transform from
+        // Nordsieck to multistep, then shifting rows to represent step advance
+        // then applying inverse transform
+        T[][] shiftedP = bigP.getData();
+        for (int i = shiftedP.length - 1; i > 0; --i) {
+            // shift rows
+            shiftedP[i] = shiftedP[i - 1];
+        }
+        shiftedP[0] = MathArrays.buildArray(field, rows);
+        Arrays.fill(shiftedP[0], field.getZero());
+        update = new Array2DRowFieldMatrix<T>(pSolver.solve(new Array2DRowFieldMatrix<T>(shiftedP, false)).getData());
+
+    }
+
+    /** Get the Nordsieck transformer for a given field and number of steps.
+     * @param field field to which the time and state vector elements belong
+     * @param nSteps number of steps of the multistep method
+     * (excluding the one being computed)
+     * @return Nordsieck transformer for the specified field and number of steps
+     * @param <T> the type of the field elements
+     */
+    public static <T extends RealFieldElement<T>> AdamsNordsieckFieldTransformer<T>
+    getInstance(final Field<T> field, final int nSteps) {
+        synchronized(CACHE) {
+            Map<Field<? extends RealFieldElement<?>>,
+                      AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>> map = CACHE.get(nSteps);
+            if (map == null) {
+                map = new HashMap<Field<? extends RealFieldElement<?>>,
+                                        AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>();
+                CACHE.put(nSteps, map);
+            }
+            @SuppressWarnings("unchecked")
+            AdamsNordsieckFieldTransformer<T> t = (AdamsNordsieckFieldTransformer<T>) map.get(field);
+            if (t == null) {
+                t = new AdamsNordsieckFieldTransformer<T>(field, nSteps);
+                map.put(field, t);
+            }
+            return t;
+
+        }
+    }
+
+    /** Build the P matrix.
+     * <p>The P matrix general terms are shifted (j+1) (-i)<sup>j</sup> terms
+     * with i being the row number starting from 1 and j being the column
+     * number starting from 1:
+     * <pre>
+     *        [  -2   3   -4    5  ... ]
+     *        [  -4  12  -32   80  ... ]
+     *   P =  [  -6  27 -108  405  ... ]
+     *        [  -8  48 -256 1280  ... ]
+     *        [          ...           ]
+     * </pre></p>
+     * @param rows number of rows of the matrix
+     * @return P matrix
+     */
+    private FieldMatrix<T> buildP(final int rows) {
+
+        final T[][] pData = MathArrays.buildArray(field, rows, rows);
+
+        for (int i = 1; i <= pData.length; ++i) {
+            // build the P matrix elements from Taylor series formulas
+            final T[] pI = pData[i - 1];
+            final int factor = -i;
+            T aj = field.getZero().add(factor);
+            for (int j = 1; j <= pI.length; ++j) {
+                pI[j - 1] = aj.multiply(j + 1);
+                aj = aj.multiply(factor);
+            }
+        }
+
+        return new Array2DRowFieldMatrix<T>(pData, false);
+
+    }
+
+    /** Initialize the high order scaled derivatives at step start.
+     * @param h step size to use for scaling
+     * @param t first steps times
+     * @param y first steps states
+     * @param yDot first steps derivatives
+     * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>,
+     * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
+     */
+
+    public Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
+                                                                   final T[][] y,
+                                                                   final T[][] yDot) {
+
+        // using Taylor series with di = ti - t0, we get:
+        //  y(ti)  - y(t0)  - di y'(t0) =   di^2 / h^2 s2 + ... +   di^k     / h^k sk + O(h^k)
+        //  y'(ti) - y'(t0)             = 2 di   / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
+        // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
+        // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
+        // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
+        // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
+        // truncated and we don't want the error terms to be included in s2 ... sk, so we need
+        // to solve also for the remainder
+        final T[][] a     = MathArrays.buildArray(field, c1.length + 1, c1.length + 1);
+        final T[][] b     = MathArrays.buildArray(field, c1.length + 1, y[0].length);
+        final T[]   y0    = y[0];
+        final T[]   yDot0 = yDot[0];
+        for (int i = 1; i < y.length; ++i) {
+
+            final T di    = t[i].subtract(t[0]);
+            final T ratio = di.divide(h);
+            T dikM1Ohk    = h.reciprocal();
+
+            // linear coefficients of equations
+            // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
+            final T[] aI    = a[2 * i - 2];
+            final T[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
+            for (int j = 0; j < aI.length; ++j) {
+                dikM1Ohk = dikM1Ohk.multiply(ratio);
+                aI[j]    = di.multiply(dikM1Ohk);
+                if (aDotI != null) {
+                    aDotI[j]  = dikM1Ohk.multiply(j + 2);
+                }
+            }
+
+            // expected value of the previous equations
+            final T[] yI    = y[i];
+            final T[] yDotI = yDot[i];
+            final T[] bI    = b[2 * i - 2];
+            final T[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
+            for (int j = 0; j < yI.length; ++j) {
+                bI[j]    = yI[j].subtract(y0[j]).subtract(di.multiply(yDot0[j]));
+                if (bDotI != null) {
+                    bDotI[j] = yDotI[j].subtract(yDot0[j]);
+                }
+            }
+
+        }
+
+        // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
+        // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
+        final FieldLUDecomposition<T> decomposition = new FieldLUDecomposition<T>(new Array2DRowFieldMatrix<T>(a, false));
+        final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<T>(b, false));
+
+        // extract just the Nordsieck vector [s2 ... sk]
+        final Array2DRowFieldMatrix<T> truncatedX =
+                        new Array2DRowFieldMatrix<T>(field, x.getRowDimension() - 1, x.getColumnDimension());
+        for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
+            for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
+                truncatedX.setEntry(i, j, x.getEntry(i, j));
+            }
+        }
+        return truncatedX;
+
+    }
+
+    /** Update the high order scaled derivatives for Adams integrators (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[], Array2DRowFieldMatrix)
+     */
+    public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
+        return update.multiply(highOrder);
+    }
+
+    /** Update the high order scaled derivatives Adams integrators (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(Array2DRowFieldMatrix)
+     */
+    public void updateHighOrderDerivativesPhase2(final T[] start,
+                                                 final T[] end,
+                                                 final Array2DRowFieldMatrix<T> highOrder) {
+        final T[][] data = highOrder.getDataRef();
+        for (int i = 0; i < data.length; ++i) {
+            final T[] dataI = data[i];
+            final T c1I = c1[i];
+            for (int j = 0; j < dataI.length; ++j) {
+                dataI[j] = dataI[j].add(c1I.multiply(start[j].subtract(end[j])));
+            }
+        }
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
new file mode 100644
index 0000000..63ab60d
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
@@ -0,0 +1,261 @@
+/*
+ * 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.math3.ode.nonstiff;
+
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.exception.MathIllegalStateException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.AbstractFieldIntegrator;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.FirstOrderFieldIntegrator;
+import org.apache.commons.math3.ode.MultistepFieldIntegrator;
+import org.apache.commons.math3.ode.TestFieldProblem1;
+import org.apache.commons.math3.ode.TestFieldProblem5;
+import org.apache.commons.math3.ode.TestFieldProblem6;
+import org.apache.commons.math3.ode.TestFieldProblemAbstract;
+import org.apache.commons.math3.ode.TestFieldProblemHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public abstract class AbstractAdamsFieldIntegratorTest {
+
+    protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
+    createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
+                     final double scalAbsoluteTolerance, final double scalRelativeTolerance);
+
+    protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
+    createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
+                     final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
+
+    @Test(expected=NumberIsTooSmallException.class)
+    public abstract void testMinStep();
+
+    protected <T extends RealFieldElement<T>> void doDimensionCheck(final Field<T> field) {
+        TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+
+        double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
+        double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+        double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
+        double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
+
+        FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
+                                                              vecAbsoluteTolerance,
+                                                              vecRelativeTolerance);
+        TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+    }
+
+    @Test
+    public abstract void testIncreasingTolerance();
+
+    protected <T extends RealFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
+                                                                             int ratioMin, int ratioMax) {
+
+        int previousCalls = Integer.MAX_VALUE;
+        for (int i = -12; i < -5; ++i) {
+            TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+            double minStep = 0;
+            double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+            double scalAbsoluteTolerance = FastMath.pow(10.0, i);
+            double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+            FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
+                                                                  scalAbsoluteTolerance,
+                                                                  scalRelativeTolerance);
+            TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+            integ.addStepHandler(handler);
+            integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+            Assert.assertTrue(handler.getMaximalValueError().getReal() > ratioMin * scalAbsoluteTolerance);
+            Assert.assertTrue(handler.getMaximalValueError().getReal() < ratioMax * scalAbsoluteTolerance);
+
+            int calls = pb.getCalls();
+            Assert.assertEquals(integ.getEvaluations(), calls);
+            Assert.assertTrue(calls <= previousCalls);
+            previousCalls = calls;
+
+        }
+
+    }
+
+    @Test(expected = MaxCountExceededException.class)
+    public abstract void exceedMaxEvaluations();
+
+    protected <T extends RealFieldElement<T>> void doExceedMaxEvaluations(final Field<T> field) {
+
+        TestFieldProblem1<T> pb  = new TestFieldProblem1<T>(field);
+        double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+
+        FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 2, 0, range, 1.0e-12, 1.0e-12);
+        TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+        integ.addStepHandler(handler);
+        integ.setMaxEvaluations(650);
+        integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+    }
+
+    @Test
+    public abstract void backward();
+
+    protected <T extends RealFieldElement<T>> void doBackward(final Field<T> field,
+                                                              final double epsilonLast,
+                                                              final double epsilonMaxValue,
+                                                              final double epsilonMaxTime,
+                                                              final String name) {
+
+        TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
+        double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+
+        AdamsFieldIntegrator<T> integ = createIntegrator(field, 4, 0, range, 1.0e-12, 1.0e-12);
+        integ.setStarterIntegrator(new PerfectStarter<T>(pb, (integ.getNSteps() + 5) / 2));
+        TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+        Assert.assertEquals(0.0, handler.getLastError().getReal(), epsilonLast);
+        Assert.assertEquals(0.0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
+        Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
+        Assert.assertEquals(name, integ.getName());
+    }
+
+    @Test
+    public abstract void polynomial();
+
+    protected <T extends RealFieldElement<T>> void doPolynomial(final Field<T> field,
+                                                                final int nLimit,
+                                                                final double epsilonBad,
+                                                                final double epsilonGood) {
+        TestFieldProblem6<T> pb = new TestFieldProblem6<T>(field);
+        double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).abs().getReal();
+
+        for (int nSteps = 2; nSteps < 8; ++nSteps) {
+            AdamsFieldIntegrator<T> integ = createIntegrator(field, nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
+            integ.setStarterIntegrator(new PerfectStarter<T>(pb, nSteps));
+            TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+            integ.addStepHandler(handler);
+            integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+            if (nSteps < nLimit) {
+                Assert.assertTrue(handler.getMaximalValueError().getReal() > epsilonBad);
+            } else {
+                Assert.assertTrue(handler.getMaximalValueError().getReal() < epsilonGood);
+            }
+        }
+
+    }
+
+    @Test(expected=MathIllegalStateException.class)
+    public abstract void testStartFailure();
+
+    protected <T extends RealFieldElement<T>> void doTestStartFailure(final Field<T> field) {
+        TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+        double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0001).getReal();
+        double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+        double scalAbsoluteTolerance = 1.0e-6;
+        double scalRelativeTolerance = 1.0e-7;
+
+        MultistepFieldIntegrator<T> integ = createIntegrator(field, 6, minStep, maxStep,
+                                                             scalAbsoluteTolerance,
+                                                             scalRelativeTolerance);
+        integ.setStarterIntegrator(new DormandPrince853FieldIntegrator<T>(field, maxStep * 0.5, maxStep, 0.1, 0.1));
+        TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+        integ.addStepHandler(handler);
+        integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+    }
+
+    private static class PerfectStarter<T extends RealFieldElement<T>> extends AbstractFieldIntegrator<T> {
+
+        private final PerfectInterpolator<T> interpolator;
+        private final int nbSteps;
+
+        public PerfectStarter(final TestFieldProblemAbstract<T> problem, final int nbSteps) {
+            super(problem.getField(), "perfect-starter");
+            this.interpolator = new PerfectInterpolator<T>(problem);
+            this.nbSteps      = nbSteps;
+        }
+
+        public FieldODEStateAndDerivative<T> integrate(FieldExpandableODE<T> equations,
+                                                       FieldODEState<T> initialState, T finalTime) {
+            T tStart = initialState.getTime().add(finalTime.subtract(initialState.getTime()).multiply(0.01));
+            getEvaluationsCounter().increment(nbSteps);
+            interpolator.setCurrentTime(initialState.getTime());
+            for (int i = 0; i < nbSteps; ++i) {
+                T tK = initialState.getTime().multiply(nbSteps - 1 - (i + 1)).add(tStart.multiply(i + 1)).divide(nbSteps - 1);
+                interpolator.setPreviousTime(interpolator.getCurrentTime());
+                interpolator.setCurrentTime(tK);
+                for (FieldStepHandler<T> handler : getStepHandlers()) {
+                    handler.handleStep(interpolator, i == nbSteps - 1);
+                }
+            }
+            return interpolator.getInterpolatedState(tStart);
+        }
+
+    }
+
+    private static class PerfectInterpolator<T extends RealFieldElement<T>> implements FieldStepInterpolator<T> {
+        private final TestFieldProblemAbstract<T> problem;
+        private T previousTime;
+        private T currentTime;
+
+        public PerfectInterpolator(final TestFieldProblemAbstract<T> problem) {
+            this.problem = problem;
+        }
+
+        public void setPreviousTime(T previousTime) {
+            this.previousTime = previousTime;
+        }
+
+        public void setCurrentTime(T currentTime) {
+            this.currentTime = currentTime;
+        }
+
+        public T getCurrentTime() {
+            return currentTime;
+        }
+
+        public boolean isForward() {
+            return problem.getFinalTime().subtract(problem.getInitialState().getTime()).getReal() >= 0;
+        }
+
+        public FieldODEStateAndDerivative<T> getPreviousState() {
+            return getInterpolatedState(previousTime);
+        }
+
+        public FieldODEStateAndDerivative<T> getCurrentState() {
+            return getInterpolatedState(currentTime);
+        }
+
+        public FieldODEStateAndDerivative<T> getInterpolatedState(T time) {
+            T[] y    = problem.computeTheoreticalState(time);
+            T[] yDot = problem.computeDerivatives(time, y);
+            return new FieldODEStateAndDerivative<T>(time, y, yDot);
+        }
+
+    }
+
+}


Mime
View raw message