commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject [3/4] [math] One step towards immutable step interpolators.
Date Tue, 01 Dec 2015 11:23:47 GMT
One step towards immutable step interpolators.

The interpolators do not expect anymore the y and yDot arrays to
be shared with integrator and be updated by it.

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

Branch: refs/heads/field-ode
Commit: 07e2b9447b5950721f39b7d027a0ca95d0cd950e
Parents: 8180c9e
Author: Luc Maisonobe <luc@apache.org>
Authored: Tue Dec 1 12:22:33 2015 +0100
Committer: Luc Maisonobe <luc@apache.org>
Committed: Tue Dec 1 12:22:33 2015 +0100

----------------------------------------------------------------------
 .../ClassicalRungeKuttaFieldIntegrator.java     |   7 +-
 ...lassicalRungeKuttaFieldStepInterpolator.java |  60 +-
 .../DormandPrince54FieldIntegrator.java         |  11 +-
 .../DormandPrince54FieldStepInterpolator.java   | 274 ++++----
 .../DormandPrince853FieldIntegrator.java        |  74 ++-
 .../DormandPrince853FieldStepInterpolator.java  | 630 +++++--------------
 .../EmbeddedRungeKuttaFieldIntegrator.java      |  31 +-
 .../ode/nonstiff/EulerFieldIntegrator.java      |   7 +-
 .../nonstiff/EulerFieldStepInterpolator.java    |  28 +-
 .../math3/ode/nonstiff/GillFieldIntegrator.java |   7 +-
 .../ode/nonstiff/GillFieldStepInterpolator.java |  62 +-
 .../nonstiff/HighamHall54FieldIntegrator.java   |  11 +-
 .../HighamHall54FieldStepInterpolator.java      |  61 +-
 .../ode/nonstiff/LutherFieldIntegrator.java     |   7 +-
 .../nonstiff/LutherFieldStepInterpolator.java   | 108 +---
 .../ode/nonstiff/MidpointFieldIntegrator.java   |   7 +-
 .../nonstiff/MidpointFieldStepInterpolator.java |  40 +-
 .../ode/nonstiff/RungeKuttaFieldIntegrator.java |  13 +-
 .../RungeKuttaFieldStepInterpolator.java        |  89 ++-
 .../nonstiff/ThreeEighthesFieldIntegrator.java  |   7 +-
 .../ThreeEighthesFieldStepInterpolator.java     |  48 +-
 .../sampling/AbstractFieldStepInterpolator.java |  21 +-
 .../commons/math3/ode/TestFieldProblem1.java    |   5 +-
 .../commons/math3/ode/TestFieldProblem5.java    |   2 +-
 .../math3/ode/TestFieldProblemAbstract.java     |  18 +-
 .../math3/ode/TestFieldProblemHandler.java      |   2 +-
 .../math3/ode/nonstiff/EulerIntegratorTest.java |   5 +-
 27 files changed, 593 insertions(+), 1042 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
index eab1bbf..064b56d 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
 
 import org.apache.commons.math3.Field;
 import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.util.MathArrays;
 
@@ -101,10 +100,8 @@ public class ClassicalRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
     /** {@inheritDoc} */
     @Override
     protected ClassicalRungeKuttaFieldStepInterpolator<T>
-        createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
-                           final T[][] yDotArray, final boolean forward,
-                           final FieldEquationsMapper<T> mapper) {
-        return new ClassicalRungeKuttaFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+        createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+        return new ClassicalRungeKuttaFieldStepInterpolator<T>(this, forward, mapper);
     }
 
 }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
index 5bea4b9..fc56870 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
 import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
 
 /**
  * This class implements a step interpolator for the classical fourth
@@ -62,17 +61,13 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
 
     /** Simple constructor.
      * @param rkIntegrator integrator being used
-     * @param y reference to the integrator array holding the state at
-     * the end of the step
-     * @param yDotArray reference to the integrator array holding all the
-     * intermediate slopes
      * @param forward integration direction indicator
      * @param mapper equations mapper for the all equations
      */
     ClassicalRungeKuttaFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
-                                             final T[] y, final T[][] yDotArray, final boolean forward,
+                                             final boolean forward,
                                              final FieldEquationsMapper<T> mapper) {
-        super(rkIntegrator, y, yDotArray, forward, mapper);
+        super(rkIntegrator, forward, mapper);
     }
 
     /** Copy constructor.
@@ -91,6 +86,7 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
     }
 
     /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
     @Override
     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
                                                                                    final T time, final T theta,
@@ -102,42 +98,28 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
         final T coeffDot1                 = oneMinusTheta.multiply(oneMinus2Theta);
         final T coeffDot23                = theta.multiply(oneMinusTheta).multiply(2);
         final T coeffDot4                 = theta.multiply(oneMinus2Theta).negate();
-        final T[] interpolatedState       = MathArrays.buildArray(theta.getField(), previousState.length);
-        final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
+        final T[] interpolatedState;
+        final T[] interpolatedDerivatives;
 
-        if ((previousState != null) && (theta.getReal() <= 0.5)) {
-            final T fourTheta2    = theta.multiply(theta).multiply(4);
-            final T s             = theta.multiply(h).divide(6.0);
-            final T coeff1        = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
-            final T coeff23       = s.multiply(theta.multiply(6).subtract(fourTheta2));
-            final T coeff4        = s.multiply(fourTheta2.subtract(theta.multiply(3)));
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                final T yDot1  = yDotK[0][i];
-                final T yDot23 = yDotK[1][i].add(yDotK[2][i]);
-                final T yDot4  = yDotK[3][i];
-                interpolatedState[i] =
-                        previousState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4));
-                interpolatedDerivatives[i] =
-                        coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4));
-            }
+        if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+            final T fourTheta2      = theta.multiply(theta).multiply(4);
+            final T s               = theta.multiply(h).divide(6.0);
+            final T coeff1          = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
+            final T coeff23         = s.multiply(theta.multiply(6).subtract(fourTheta2));
+            final T coeff4          = s.multiply(fourTheta2.subtract(theta.multiply(3)));
+            interpolatedState       = previousStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
+            interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
         } else {
-            final T fourTheta     = theta.multiply(4);
-            final T s             = oneMinusThetaH.divide(6);
-            final T coeff1        = s.multiply(theta.multiply(fourTheta.negate().add(5)).subtract(1));
-            final T coeff23       = s.multiply(theta.multiply(fourTheta.subtract(2)).subtract(2));
-            final T coeff4        = s.multiply(theta.multiply(fourTheta.negate().subtract(1)).subtract(1));
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                final T yDot1  = yDotK[0][i];
-                final T yDot23 = yDotK[1][i].add(yDotK[2][i]);
-                final T yDot4  = yDotK[3][i];
-                interpolatedState[i] =
-                        currentState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4));
-                interpolatedDerivatives[i] =
-                        coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4));
-            }
+            final T fourTheta       = theta.multiply(4);
+            final T s               = oneMinusThetaH.divide(6);
+            final T coeff1          = s.multiply(theta.multiply(fourTheta.negate().add(5)).subtract(1));
+            final T coeff23         = s.multiply(theta.multiply(fourTheta.subtract(2)).subtract(2));
+            final T coeff4          = s.multiply(theta.multiply(fourTheta.negate().subtract(1)).subtract(1));
+            interpolatedState       = currentStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
+            interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
         }
 
-        return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+        return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
index 7475c66..09bbc4a 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
 
 import org.apache.commons.math3.Field;
 import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.util.MathArrays;
 import org.apache.commons.math3.util.MathUtils;
@@ -93,7 +92,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
                                           final double minStep, final double maxStep,
                                           final double scalAbsoluteTolerance,
                                           final double scalRelativeTolerance) {
-        super(field, METHOD_NAME, true,
+        super(field, METHOD_NAME, 6,
               minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
         e1 = fraction(    71,  57600);
         e3 = fraction(   -71,  16695);
@@ -119,7 +118,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
                                           final double minStep, final double maxStep,
                                           final double[] vecAbsoluteTolerance,
                                           final double[] vecRelativeTolerance) {
-        super(field, METHOD_NAME, true,
+        super(field, METHOD_NAME, 6,
               minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
         e1 = fraction(    71,  57600);
         e3 = fraction(   -71,  16695);
@@ -190,10 +189,8 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
     /** {@inheritDoc} */
     @Override
     protected DormandPrince54FieldStepInterpolator<T>
-        createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
-                           final T[][] yDotArray, final boolean forward,
-                           final FieldEquationsMapper<T> mapper) {
-        return new DormandPrince54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+        createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+        return new DormandPrince54FieldStepInterpolator<T>(this, forward, mapper);
     }
 
     /** {@inheritDoc} */

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
index d3fb208..23ba0e9 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
 import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
 
 /**
  * This class represents an interpolator over the last step during an
@@ -37,75 +36,63 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
       extends RungeKuttaFieldStepInterpolator<T> {
 
     /** Last row of the Butcher-array internal weights, element 0. */
-    private static final double A70 =    35.0 /  384.0;
+    private final T a70;
 
     // element 1 is zero, so it is neither stored nor used
 
     /** Last row of the Butcher-array internal weights, element 2. */
-    private static final double A72 =   500.0 / 1113.0;
+    private final T a72;
 
     /** Last row of the Butcher-array internal weights, element 3. */
-    private static final double A73 =   125.0 /  192.0;
+    private final T a73;
 
     /** Last row of the Butcher-array internal weights, element 4. */
-    private static final double A74 = -2187.0 / 6784.0;
+    private final T a74;
 
     /** Last row of the Butcher-array internal weights, element 5. */
-    private static final double A75 =    11.0 /   84.0;
+    private final T a75;
 
     /** Shampine (1986) Dense output, element 0. */
-    private static final double D0 =  -12715105075.0 /  11282082432.0;
+    private final T d0;
 
     // element 1 is zero, so it is neither stored nor used
 
     /** Shampine (1986) Dense output, element 2. */
-    private static final double D2 =   87487479700.0 /  32700410799.0;
+    private final T d2;
 
     /** Shampine (1986) Dense output, element 3. */
-    private static final double D3 =  -10690763975.0 /   1880347072.0;
+    private final T d3;
 
     /** Shampine (1986) Dense output, element 4. */
-    private static final double D4 =  701980252875.0 / 199316789632.0;
+    private final T d4;
 
     /** Shampine (1986) Dense output, element 5. */
-    private static final double D5 =   -1453857185.0 /    822651844.0;
+    private final T d5;
 
     /** Shampine (1986) Dense output, element 6. */
-    private static final double D6 =      69997945.0 /     29380423.0;
-
-    /** First vector for interpolation. */
-    private T[] v1;
-
-    /** Second vector for interpolation. */
-    private T[] v2;
-
-    /** Third vector for interpolation. */
-    private T[] v3;
-
-    /** Fourth vector for interpolation. */
-    private T[] v4;
-
-    /** Initialization indicator for the interpolation vectors. */
-    private boolean vectorsInitialized;
+    private final T d6;
 
     /** Simple constructor.
      * @param rkIntegrator integrator being used
-     * @param y reference to the integrator array holding the state at
-     * the end of the step
-     * @param yDotArray reference to the integrator array holding all the
-     * intermediate slopes
      * @param forward integration direction indicator
      * @param mapper equations mapper for the all equations
      */
     DormandPrince54FieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
-                                         final T[] y, final T[][] yDotArray, final boolean forward,
+                                         final boolean forward,
                                          final FieldEquationsMapper<T> mapper) {
-        super(rkIntegrator, y, yDotArray, forward, mapper);
-        v1 = null;
-        v2 = null;
-        v3 = null;
-        v4 = null;
-        vectorsInitialized = false;
+        super(rkIntegrator, forward, mapper);
+        final T one = rkIntegrator.getField().getOne();
+        a70 = one.multiply(   35.0).divide( 384.0);
+        a72 = one.multiply(  500.0).divide(1113.0);
+        a73 = one.multiply(  125.0).divide( 192.0);
+        a74 = one.multiply(-2187.0).divide(6784.0);
+        a75 = one.multiply(   11.0).divide(  84.0);
+        d0  = one.multiply(-12715105075.0).divide( 11282082432.0);
+        d2  = one.multiply( 87487479700.0).divide( 32700410799.0);
+        d3  = one.multiply(-10690763975.0).divide(  1880347072.0);
+        d4  = one.multiply(701980252875.0).divide(199316789632.0);
+        d5  = one.multiply( -1453857185.0).divide(   822651844.0);
+        d6  = one.multiply(    69997945.0).divide(    29380423.0);
     }
 
     /** Copy constructor.
@@ -116,24 +103,17 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
     DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T> interpolator) {
 
         super(interpolator);
-
-        if (interpolator.v1 == null) {
-
-            v1 = null;
-            v2 = null;
-            v3 = null;
-            v4 = null;
-            vectorsInitialized = false;
-
-        } else {
-
-            v1 = interpolator.v1.clone();
-            v2 = interpolator.v2.clone();
-            v3 = interpolator.v3.clone();
-            v4 = interpolator.v4.clone();
-            vectorsInitialized = interpolator.vectorsInitialized;
-
-        }
+        a70 = interpolator.a70;
+        a72 = interpolator.a72;
+        a73 = interpolator.a73;
+        a74 = interpolator.a74;
+        a75 = interpolator.a75;
+        d0  = interpolator.d0;
+        d2  = interpolator.d2;
+        d3  = interpolator.d3;
+        d4  = interpolator.d4;
+        d5  = interpolator.d5;
+        d6  = interpolator.d6;
 
     }
 
@@ -143,58 +123,13 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
         return new DormandPrince54FieldStepInterpolator<T>(this);
     }
 
-
-    /** {@inheritDoc} */
-    @Override
-    public void storeState(final FieldODEStateAndDerivative<T> state) {
-        super.storeState(state);
-        vectorsInitialized = false;
-    }
-
     /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
     @Override
     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
                                                                                    final T time, final T theta,
                                                                                    final T oneMinusThetaH) {
 
-        if (! vectorsInitialized) {
-
-            if (v1 == null) {
-                v1 = MathArrays.buildArray(time.getField(), previousState.length);
-                v2 = MathArrays.buildArray(time.getField(), previousState.length);
-                v3 = MathArrays.buildArray(time.getField(), previousState.length);
-                v4 = MathArrays.buildArray(time.getField(), previousState.length);
-            }
-
-            // no step finalization is needed for this interpolator
-
-            // we need to compute the interpolation vectors for this time step
-            for (int i = 0; i < previousState.length; ++i) {
-                final T yDot0 = yDotK[0][i];
-                final T yDot2 = yDotK[2][i];
-                final T yDot3 = yDotK[3][i];
-                final T yDot4 = yDotK[4][i];
-                final T yDot5 = yDotK[5][i];
-                final T yDot6 = yDotK[6][i];
-                v1[i] =     yDot0.multiply(A70).
-                        add(yDot2.multiply(A72)).
-                        add(yDot3.multiply(A73)).
-                        add(yDot4.multiply(A74)).
-                        add(yDot5.multiply(A75));
-                v2[i] = yDot0.subtract(v1[i]);
-                v3[i] = v1[i].subtract(v2[i]).subtract(yDot6);
-                v4[i] =     yDot0.multiply(D0).
-                        add(yDot2.multiply(D2)).
-                        add(yDot3.multiply(D3)).
-                        add(yDot4.multiply(D4)).
-                        add(yDot5.multiply(D5)).
-                        add(yDot6.multiply(D6));
-            }
-
-            vectorsInitialized = true;
-
-        }
-
         // interpolate
         final T one      = theta.getField().getOne();
         final T eta      = one.subtract(theta);
@@ -202,35 +137,116 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
         final T dot2     = one.subtract(twoTheta);
         final T dot3     = theta.multiply(theta.multiply(-3).add(2));
         final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
-        final T[] interpolatedState       = MathArrays.buildArray(theta.getField(), previousState.length);
-        final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
-        if ((previousState != null) && (theta.getReal() <= 0.5)) {
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                interpolatedState[i] = previousState[i].
-                                       add(theta.multiply(h.multiply(v1[i].
-                                                                     add(eta.multiply(v2[i].
-                                                                                      add(theta.multiply(v3[i].
-                                                                                                         add(eta.multiply(v4[i])))))))));
-                interpolatedDerivatives[i] =                   v1[i].
-                                             add(dot2.multiply(v2[i])).
-                                             add(dot3.multiply(v3[i])).
-                                             add(dot4.multiply(v4[i]));
-            }
+        final T[] interpolatedState;
+        final T[] interpolatedDerivatives;
+        if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+            final T f1        = h.multiply(theta);
+            final T f2        = f1.multiply(eta);
+            final T f3        = f2.multiply(theta);
+            final T f4        = f3.multiply(eta);
+            final T coeff0    = f1.multiply(a70).
+                                subtract(f2.multiply(a70.subtract(1))).
+                                add(f3.multiply(a70.multiply(2).subtract(1))).
+                                add(f4.multiply(d0));
+            final T coeff1    = theta.getField().getZero();
+            final T coeff2    = f1.multiply(a72).
+                                subtract(f2.multiply(a72)).
+                                add(f3.multiply(a72.multiply(2))).
+                                add(f4.multiply(d2));
+            final T coeff3    = f1.multiply(a73).
+                                subtract(f2.multiply(a73)).
+                                add(f3.multiply(a73.multiply(2))).
+                                add(f4.multiply(d3));
+            final T coeff4    = f1.multiply(a74).
+                                subtract(f2.multiply(a74)).
+                                add(f3.multiply(a74.multiply(2))).
+                                add(f4.multiply(d4));
+            final T coeff5    = f1.multiply(a75).
+                                subtract(f2.multiply(a75)).
+                                add(f3.multiply(a75.multiply(2))).
+                                add(f4.multiply(d5));
+            final T coeff6    = f4.multiply(d6).subtract(f3);
+            final T coeffDot0 = a70.
+                                subtract(dot2.multiply(a70.subtract(1))).
+                                add(dot3.multiply(a70.multiply(2).subtract(1))).
+                                add(dot4.multiply(d0));
+            final T coeffDot1 = theta.getField().getZero();
+            final T coeffDot2 = a72.
+                                subtract(dot2.multiply(a72)).
+                                add(dot3.multiply(a72.multiply(2))).
+                                add(dot4.multiply(d2));
+            final T coeffDot3 = a73.
+                                subtract(dot2.multiply(a73)).
+                                add(dot3.multiply(a73.multiply(2))).
+                                add(dot4.multiply(d3));
+            final T coeffDot4 = a74.
+                                subtract(dot2.multiply(a74)).
+                                add(dot3.multiply(a74.multiply(2))).
+                                add(dot4.multiply(d4));
+            final T coeffDot5 = a75.
+                                subtract(dot2.multiply(a75)).
+                                add(dot3.multiply(a75.multiply(2))).
+                                add(dot4.multiply(d5));
+            final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
+            interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
+                                                                     coeff4, coeff5, coeff6);
+            interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
+                                                                  coeffDot4, coeffDot5, coeffDot6);
         } else {
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                interpolatedState[i] = currentState[i].
-                                       subtract(oneMinusThetaH.multiply(v1[i].
-                                                                        subtract(theta.multiply(v2[i].
-                                                                                                add(theta.multiply(v3[i].
-                                                                                                                   add(eta.multiply(v4[i]))))))));
-                interpolatedDerivatives[i] =                   v1[i].
-                                             add(dot2.multiply(v2[i])).
-                                             add(dot3.multiply(v3[i])).
-                                             add(dot4.multiply(v4[i]));
-            }
+            final T f1        = oneMinusThetaH.negate();
+            final T f2        = oneMinusThetaH.multiply(theta);
+            final T f3        = f2.multiply(theta);
+            final T f4        = f3.multiply(eta);
+            final T coeff0    = f1.multiply(a70).
+                                subtract(f2.multiply(a70.subtract(1))).
+                                add(f3.multiply(a70.multiply(2).subtract(1))).
+                                add(f4.multiply(d0));
+            final T coeff1    = theta.getField().getZero();
+            final T coeff2    = f1.multiply(a72).
+                                subtract(f2.multiply(a72)).
+                                add(f3.multiply(a72.multiply(2))).
+                                add(f4.multiply(d2));
+            final T coeff3    = f1.multiply(a73).
+                                subtract(f2.multiply(a73)).
+                                add(f3.multiply(a73.multiply(2))).
+                                add(f4.multiply(d3));
+            final T coeff4    = f1.multiply(a74).
+                                subtract(f2.multiply(a74)).
+                                add(f3.multiply(a74.multiply(2))).
+                                add(f4.multiply(d4));
+            final T coeff5    = f1.multiply(a75).
+                                subtract(f2.multiply(a75)).
+                                add(f3.multiply(a75.multiply(2))).
+                                add(f4.multiply(d5));
+            final T coeff6    = f4.multiply(d6).subtract(f3);
+            final T coeffDot0 = a70.
+                                subtract(dot2.multiply(a70.subtract(1))).
+                                add(dot3.multiply(a70.multiply(2).subtract(1))).
+                                add(dot4.multiply(d0));
+            final T coeffDot1 = theta.getField().getZero();
+            final T coeffDot2 = a72.
+                                subtract(dot2.multiply(a72)).
+                                add(dot3.multiply(a72.multiply(2))).
+                                add(dot4.multiply(d2));
+            final T coeffDot3 = a73.
+                                subtract(dot2.multiply(a73)).
+                                add(dot3.multiply(a73.multiply(2))).
+                                add(dot4.multiply(d3));
+            final T coeffDot4 = a74.
+                                subtract(dot2.multiply(a74)).
+                                add(dot3.multiply(a74.multiply(2))).
+                                add(dot4.multiply(d4));
+            final T coeffDot5 = a75.
+                                subtract(dot2.multiply(a75)).
+                                add(dot3.multiply(a75.multiply(2))).
+                                add(dot4.multiply(d5));
+            final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
+            interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
+                                                                    coeff4, coeff5, coeff6);
+            interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
+                                                                  coeffDot4, coeffDot5, coeffDot6);
         }
-
-        return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+        return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
index 6f6b436..569cba4 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
 
 import org.apache.commons.math3.Field;
 import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.util.MathArrays;
 import org.apache.commons.math3.util.MathUtils;
@@ -134,7 +133,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
                                            final double minStep, final double maxStep,
                                            final double scalAbsoluteTolerance,
                                            final double scalRelativeTolerance) {
-        super(field, METHOD_NAME, true,
+        super(field, METHOD_NAME, 12,
               minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
         e1_01 = fraction(        116092271.0,       8848465920.0);
         e1_06 = fraction(         -1871647.0,          1527680.0);
@@ -170,7 +169,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
                                            final double minStep, final double maxStep,
                                            final double[] vecAbsoluteTolerance,
                                            final double[] vecRelativeTolerance) {
-        super(field, METHOD_NAME, true,
+        super(field, METHOD_NAME, 12,
               minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
         e1_01 = fraction(        116092271.0,       8848465920.0);
         e1_06 = fraction(         -1871647.0,          1527680.0);
@@ -196,7 +195,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
 
         final T sqrt6 = getField().getOne().multiply(6).sqrt();
 
-        final T[] c = MathArrays.buildArray(getField(), 12);
+        final T[] c = MathArrays.buildArray(getField(), 15);
         c[ 0] = sqrt6.add(-6).divide(-67.5);
         c[ 1] = sqrt6.add(-6).divide(-45.0);
         c[ 2] = sqrt6.add(-6).divide(-30.0);
@@ -209,6 +208,9 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
         c[ 9] = fraction(6, 7);
         c[10] = getField().getOne();
         c[11] = getField().getOne();
+        c[12] = fraction(1.0, 10.0);
+        c[13] = fraction(1.0, 5.0);
+        c[14] = fraction(7.0, 9.0);
 
         return c;
 
@@ -220,7 +222,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
 
         final T sqrt6 = getField().getOne().multiply(6).sqrt();
 
-        final T[][] a = MathArrays.buildArray(getField(), 12, -1);
+        final T[][] a = MathArrays.buildArray(getField(), 15, -1);
         for (int i = 0; i < a.length; ++i) {
             a[i] = MathArrays.buildArray(getField(), i + 1);
         }
@@ -302,9 +304,8 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
         a[10][ 9] = fraction(226716250.0, 18341897.0);
         a[10][10] = fraction(1371316744.0, 2131383595.0);
 
-        // this stage should be for interpolation only, but since it is the same
-        // stage as the first evaluation of the next step, we perform it
-        // here at no cost by specifying this is an fsal method
+        // the following stage is both for interpolation and the first stage in next step
+        // (the coefficients are identical to the B array)
         a[11][ 0] = fraction(104257.0, 1920240.0);
         a[11][ 1] = getField().getZero();
         a[11][ 2] = getField().getZero();
@@ -318,6 +319,52 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
         a[11][10] = fraction(171414593.0, 851261400.0);
         a[11][11] = fraction(137909.0, 3084480.0);
 
+        // the following stages are for interpolation only
+        a[12][ 0] = fraction(      13481885573.0, 240030000000.0)     .subtract(a[11][0]);
+        a[12][ 1] = getField().getZero();
+        a[12][ 2] = getField().getZero();
+        a[12][ 3] = getField().getZero();
+        a[12][ 4] = getField().getZero();
+        a[12][ 5] = getField().getZero()                              .subtract(a[11][5]);
+        a[12][ 6] = fraction(     139418837528.0, 549975234375.0)     .subtract(a[11][6]);
+        a[12][ 7] = fraction(  -11108320068443.0, 45111937500000.0)   .subtract(a[11][7]);
+        a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0).subtract(a[11][8]);
+        a[12][ 9] = fraction(         57799439.0, 377055000.0)        .subtract(a[11][9]);
+        a[12][10] = fraction(     793322643029.0, 96734250000000.0)   .subtract(a[11][10]);
+        a[12][11] = fraction(       1458939311.0, 192780000000.0)     .subtract(a[11][11]);
+        a[12][12]  = fraction(            -4149.0, 500000.0);
+
+        a[13][ 0] = fraction(    1595561272731.0, 50120273500000.0)   .subtract(a[11][0]);
+        a[13][ 1] = getField().getZero();
+        a[13][ 2] = getField().getZero();
+        a[13][ 3] = getField().getZero();
+        a[13][ 4] = getField().getZero();
+        a[13][ 5] = fraction(     975183916491.0, 34457688031250.0)   .subtract(a[11][5]);
+        a[13][ 6] = fraction(   38492013932672.0, 718912673015625.0)  .subtract(a[11][6]);
+        a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0).subtract(a[11][7]);
+        a[13][ 8] = getField().getZero()                              .subtract(a[11][8]);
+        a[13][ 9] = getField().getZero()                              .subtract(a[11][9]);
+        a[13][10] = fraction(   -2538710946863.0, 23431227861250000.0).subtract(a[11][10]);
+        a[13][11] = fraction(       8824659001.0, 23066716781250.0)   .subtract(a[11][11]);
+        a[13][12] = fraction(     -11518334563.0, 33831184612500.0);
+        a[13][13] = fraction(       1912306948.0, 13532473845.0);
+
+        a[14][ 0] = fraction(     -13613986967.0, 31741908048.0)      .subtract(a[11][0]);
+        a[14][ 1] = getField().getZero();
+        a[14][ 2] = getField().getZero();
+        a[14][ 3] = getField().getZero();
+        a[14][ 4] = getField().getZero();
+        a[14][ 5] = fraction(      -4755612631.0, 1012344804.0)       .subtract(a[11][5]);
+        a[14][ 6] = fraction(   42939257944576.0, 5588559685701.0)    .subtract(a[11][6]);
+        a[14][ 7] = fraction(   77881972900277.0, 19140370552944.0)   .subtract(a[11][7]);
+        a[14][ 8] = fraction(   22719829234375.0, 63689648654052.0)   .subtract(a[11][8]);
+        a[14][ 9] = getField().getZero()                              .subtract(a[11][9]);
+        a[14][10] = getField().getZero()                              .subtract(a[11][10]);
+        a[14][11] = getField().getZero()                              .subtract(a[11][11]);
+        a[14][12] = fraction(      -1199007803.0, 857031517296.0);
+        a[14][13] = fraction(     157882067000.0, 53564469831.0);
+        a[14][14] = fraction(    -290468882375.0, 31741908048.0);
+
         return a;
 
     }
@@ -325,7 +372,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
     /** {@inheritDoc} */
     @Override
     protected T[] getB() {
-        final T[] b = MathArrays.buildArray(getField(), 13);
+        final T[] b = MathArrays.buildArray(getField(), 16);
         b[ 0] = fraction(104257, 1929240);
         b[ 1] = getField().getZero();
         b[ 2] = getField().getZero();
@@ -339,16 +386,17 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
         b[10] = fraction(      171414593.0,       851261400.0);
         b[11] = fraction(         137909.0,         3084480.0);
         b[12] = getField().getZero();
+        b[13] = getField().getZero();
+        b[14] = getField().getZero();
+        b[15] = getField().getZero();
         return b;
     }
 
     /** {@inheritDoc} */
     @Override
     protected DormandPrince853FieldStepInterpolator<T>
-        createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
-                           final T[][] yDotArray, final boolean forward,
-                           final FieldEquationsMapper<T> mapper) {
-        return new DormandPrince853FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+        createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+        return new DormandPrince853FieldStepInterpolator<T>(this, forward, mapper);
     }
 
     /** {@inheritDoc} */

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
index 3966121..7c1cdf5 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
@@ -17,7 +17,6 @@
 
 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.MaxCountExceededException;
 import org.apache.commons.math3.ode.AbstractFieldIntegrator;
@@ -38,268 +37,143 @@ import org.apache.commons.math3.util.MathArrays;
 class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
     extends RungeKuttaFieldStepInterpolator<T> {
 
-    /** Propagation weights, element 1. */
-    private final T b_01;
-
-    // elements 2 to 5 are zero, so they are neither stored nor used
-
-    /** Propagation weights, element 6. */
-    private final T b_06;
-
-    /** Propagation weights, element 7. */
-    private final T b_07;
-
-    /** Propagation weights, element 8. */
-    private final T b_08;
-
-    /** Propagation weights, element 9. */
-    private final T b_09;
-
-    /** Propagation weights, element 10. */
-    private final T b_10;
-
-    /** Propagation weights, element 11. */
-    private final T b_11;
-
-    /** Propagation weights, element 12. */
-    private final T b_12;
-
-    /** Time step for stage 14 (interpolation only). */
-    private final T c14;
-
-    /** Internal weights for stage 14, element 1. */
-    private final T k14_01;
-
-    // elements 2 to 5 are zero, so they are neither stored nor used
-
-    /** Internal weights for stage 14, element 6. */
-    private final T k14_06;
-
-    /** Internal weights for stage 14, element 7. */
-    private final T k14_07;
-
-    /** Internal weights for stage 14, element 8. */
-    private final T k14_08;
-
-    /** Internal weights for stage 14, element 9. */
-    private final T k14_09;
-
-    /** Internal weights for stage 14, element 10. */
-    private final T k14_10;
-
-    /** Internal weights for stage 14, element 11. */
-    private final T k14_11;
-
-    /** Internal weights for stage 14, element 12. */
-    private final T k14_12;
-
-    /** Internal weights for stage 14, element 13. */
-    private final T k14_13;
-
-    /** Time step for stage 15 (interpolation only). */
-    private final T c15;
-
-
-    /** Internal weights for stage 15, element 1. */
-    private final T k15_01;
-
-    // elements 2 to 5 are zero, so they are neither stored nor used
-
-    /** Internal weights for stage 15, element 6. */
-    private final T k15_06;
-
-    /** Internal weights for stage 15, element 7. */
-    private final T k15_07;
-
-    /** Internal weights for stage 15, element 8. */
-    private final T k15_08;
-
-    /** Internal weights for stage 15, element 9. */
-    private final T k15_09;
-
-    /** Internal weights for stage 15, element 10. */
-    private final T k15_10;
-
-    /** Internal weights for stage 15, element 11. */
-    private final T k15_11;
-
-    /** Internal weights for stage 15, element 12. */
-    private final T k15_12;
-
-    /** Internal weights for stage 15, element 13. */
-    private final T k15_13;
-
-    /** Internal weights for stage 15, element 14. */
-    private final T k15_14;
-
-    /** Time step for stage 16 (interpolation only). */
-    private final T c16;
-
-
-    /** Internal weights for stage 16, element 1. */
-    private final T k16_01;
-
-    // elements 2 to 5 are zero, so they are neither stored nor used
-
-    /** Internal weights for stage 16, element 6. */
-    private final T k16_06;
-
-    /** Internal weights for stage 16, element 7. */
-    private final T k16_07;
-
-    /** Internal weights for stage 16, element 8. */
-    private final T k16_08;
-
-    /** Internal weights for stage 16, element 9. */
-    private final T k16_09;
-
-    /** Internal weights for stage 16, element 10. */
-    private final T k16_10;
-
-    /** Internal weights for stage 16, element 11. */
-    private final T k16_11;
-
-    /** Internal weights for stage 16, element 12. */
-    private final T k16_12;
-
-    /** Internal weights for stage 16, element 13. */
-    private final T k16_13;
-
-    /** Internal weights for stage 16, element 14. */
-    private final T k16_14;
-
-    /** Internal weights for stage 16, element 15. */
-    private final T k16_15;
-
     /** Interpolation weights.
      * (beware that only the non-null values are in the table)
      */
     private final T[][] d;
 
-    /** Last evaluations. */
-    private T[][] yDotKLast;
-
-    /** Vectors for interpolation. */
-    private T[][] v;
-
-    /** Initialization indicator for the interpolation vectors. */
-    private boolean vectorsInitialized;
-
     /** Simple constructor.
      * @param rkIntegrator integrator being used
-     * @param y reference to the integrator array holding the state at
-     * the end of the step
-     * @param yDotArray reference to the integrator array holding all the
-     * intermediate slopes
      * @param forward integration direction indicator
      * @param mapper equations mapper for the all equations
      */
     DormandPrince853FieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
-                                          final T[] y, final T[][] yDotArray, final boolean forward,
+                                          final boolean forward,
                                           final FieldEquationsMapper<T> mapper) {
-        super(rkIntegrator, y, yDotArray, forward, mapper);
-        yDotKLast = null;
-        v         = null;
-        vectorsInitialized = false;
-
-        b_01   = fraction(        104257.0, 1920240.0);
-        b_06   = fraction(       3399327.0, 763840.0);
-        b_07   = fraction(      66578432.0, 35198415.0);
-        b_08   = fraction(   -1674902723.0, 288716400.0);
-        b_09   = fraction(54980371265625.0, 176692375811392.0);
-        b_10   = fraction(       -734375.0, 4826304.0);
-        b_11   = fraction(     171414593.0, 851261400.0);
-        b_12   = fraction(        137909.0, 3084480.0);
-        c14    = fraction(1.0, 10.0);
-        k14_01 = fraction(      13481885573.0, 240030000000.0)     .subtract(b_01);
-        k14_06 = integrator.getField().getZero()                   .subtract(b_06);
-        k14_07 = fraction(     139418837528.0, 549975234375.0)     .subtract(b_07);
-        k14_08 = fraction(  -11108320068443.0, 45111937500000.0)   .subtract(b_08);
-        k14_09 = fraction(-1769651421925959.0, 14249385146080000.0).subtract(b_09);
-        k14_10 = fraction(         57799439.0, 377055000.0)        .subtract(b_10);
-        k14_11 = fraction(     793322643029.0, 96734250000000.0)   .subtract(b_11);
-        k14_12 = fraction(       1458939311.0, 192780000000.0)     .subtract(b_12);
-        k14_13 = fraction(            -4149.0, 500000.0);
-        c15    = fraction(1.0, 5.0);
-        k15_01 = fraction(    1595561272731.0, 50120273500000.0)   .subtract(b_01);
-        k15_06 = fraction(     975183916491.0, 34457688031250.0)   .subtract(b_06);
-        k15_07 = fraction(   38492013932672.0, 718912673015625.0)  .subtract(b_07);
-        k15_08 = fraction(-1114881286517557.0, 20298710767500000.0).subtract(b_08);
-        k15_09 = integrator.getField().getZero()                   .subtract(b_09);
-        k15_10 = integrator.getField().getZero()                   .subtract(b_10);
-        k15_11 = fraction(   -2538710946863.0, 23431227861250000.0).subtract(b_11);
-        k15_12 = fraction(       8824659001.0, 23066716781250.0)   .subtract(b_12);
-        k15_13 = fraction(     -11518334563.0, 33831184612500.0);
-        k15_14 = fraction(       1912306948.0, 13532473845.0);
-        c16    = fraction(7.0, 9.0);
-        k16_01 = fraction(     -13613986967.0, 31741908048.0)      .subtract(b_01);
-        k16_06 = fraction(      -4755612631.0, 1012344804.0)       .subtract(b_06);
-        k16_07 = fraction(   42939257944576.0, 5588559685701.0)    .subtract(b_07);
-        k16_08 = fraction(   77881972900277.0, 19140370552944.0)   .subtract(b_08);
-        k16_09 = fraction(   22719829234375.0, 63689648654052.0)   .subtract(b_09);
-        k16_10 = integrator.getField().getZero()                   .subtract(b_10);
-        k16_11 = integrator.getField().getZero()                   .subtract(b_11);
-        k16_12 = integrator.getField().getZero()                   .subtract(b_12);
-        k16_13 = fraction(      -1199007803.0, 857031517296.0);
-        k16_14 = fraction(     157882067000.0, 53564469831.0);
-        k16_15 = fraction(    -290468882375.0, 31741908048.0);
-
-        /** Interpolation weights.
-         * (beware that only the non-null values are in the table)
-         */
-        d = MathArrays.buildArray(integrator.getField(), 4, 12);
-
-         d[0][ 0] = fraction(        -17751989329.0, 2106076560.0);
-         d[0][ 1] = fraction(          4272954039.0, 7539864640.0);
-         d[0][ 2] = fraction(       -118476319744.0, 38604839385.0);
-         d[0][ 3] = fraction(        755123450731.0, 316657731600.0);
-         d[0][ 4] = fraction( 3692384461234828125.0, 1744130441634250432.0);
-         d[0][ 5] = fraction(         -4612609375.0, 5293382976.0);
-         d[0][ 6] = fraction(       2091772278379.0, 933644586600.0);
-         d[0][ 7] = fraction(          2136624137.0, 3382989120.0);
-         d[0][ 8] = fraction(             -126493.0, 1421424.0);
-         d[0][ 9] = fraction(            98350000.0, 5419179.0);
-         d[0][10] = fraction(           -18878125.0, 2053168.0);
-         d[0][11] = fraction(         -1944542619.0, 438351368.0);
-
-         d[1][ 0] = fraction(         32941697297.0, 3159114840.0);
-         d[1][ 1] = fraction(        456696183123.0, 1884966160.0);
-         d[1][ 2] = fraction(      19132610714624.0, 115814518155.0);
-         d[1][ 3] = fraction(    -177904688592943.0, 474986597400.0);
-         d[1][ 4] = fraction(-4821139941836765625.0, 218016305204281304.0);
-         d[1][ 5] = fraction(         30702015625.0, 3970037232.0);
-         d[1][ 6] = fraction(     -85916079474274.0, 2800933759800.0);
-         d[1][ 7] = fraction(         -5919468007.0, 634310460.0);
-         d[1][ 8] = fraction(             2479159.0, 157936.0);
-         d[1][ 9] = fraction(           -18750000.0, 602131.0);
-         d[1][10] = fraction(           -19203125.0, 2053168.0);
-         d[1][11] = fraction(         15700361463.0, 438351368.0);
-
-         d[2][ 0] = fraction(         12627015655.0, 631822968.0);
-         d[2][ 1] = fraction(        -72955222965.0, 188496616.0);
-         d[2][ 2] = fraction(     -13145744952320.0, 69488710893.0);
-         d[2][ 3] = fraction(      30084216194513.0, 56998391688.0);
-         d[2][ 4] = fraction( -296858761006640625.0, 25648977082856624.0);
-         d[2][ 5] = fraction(           569140625.0, 82709109.0);
-         d[2][ 6] = fraction(        -18684190637.0, 18672891732.0);
-         d[2][ 7] = fraction(            69644045.0, 89549712.0);
-         d[2][ 8] = fraction(           -11847025.0, 4264272.0);
-         d[2][ 9] = fraction(          -978650000.0, 16257537.0);
-         d[2][10] = fraction(           519371875.0, 6159504.0);
-         d[2][11] = fraction(          5256837225.0, 438351368.0);
-
-         d[3][ 0] = fraction(          -450944925.0, 17550638.0);
-         d[3][ 1] = fraction(        -14532122925.0, 94248308.0);
-         d[3][ 2] = fraction(       -595876966400.0, 2573655959.0);
-         d[3][ 3] = fraction(        188748653015.0, 527762886.0);
-         d[3][ 4] = fraction( 2545485458115234375.0, 27252038150535163.0);
-         d[3][ 5] = fraction(         -1376953125.0, 36759604.0);
-         d[3][ 6] = fraction(         53995596795.0, 518691437.0);
-         d[3][ 7] = fraction(           210311225.0, 7047894.0);
-         d[3][ 8] = fraction(            -1718875.0, 39484.0);
-         d[3][ 9] = fraction(            58000000.0, 602131.0);
-         d[3][10] = fraction(            -1546875.0, 39484.0);
-         d[3][11] = fraction(         -1262172375.0, 8429834.0);
+        super(rkIntegrator, forward, mapper);
+
+        // interpolation weights
+        d = MathArrays.buildArray(integrator.getField(), 4, 16);
+
+        // this row is the same as the b array
+        d[0][ 0] = fraction(104257, 1929240);
+        d[0][ 1] = integrator.getField().getZero();
+        d[0][ 2] = integrator.getField().getZero();
+        d[0][ 3] = integrator.getField().getZero();
+        d[0][ 4] = integrator.getField().getZero();
+        d[0][ 5] = fraction(        3399327.0,          763840.0);
+        d[0][ 6] = fraction(       66578432.0,        35198415.0);
+        d[0][ 7] = fraction(    -1674902723.0,       288716400.0);
+        d[0][ 8] = fraction( 54980371265625.0, 176692375811392.0);
+        d[0][ 9] = fraction(        -734375.0,         4826304.0);
+        d[0][10] = fraction(      171414593.0,       851261400.0);
+        d[0][11] = fraction(         137909.0,         3084480.0);
+        d[0][12] = integrator.getField().getZero();
+        d[0][13] = integrator.getField().getZero();
+        d[0][14] = integrator.getField().getZero();
+        d[0][15] = integrator.getField().getZero();
+
+        d[1][ 0] = d[0][ 0].negate().add(1);
+        d[1][ 1] = d[0][ 1].negate();
+        d[1][ 2] = d[0][ 2].negate();
+        d[1][ 3] = d[0][ 3].negate();
+        d[1][ 4] = d[0][ 4].negate();
+        d[1][ 5] = d[0][ 5].negate();
+        d[1][ 6] = d[0][ 6].negate();
+        d[1][ 7] = d[0][ 7].negate();
+        d[1][ 8] = d[0][ 8].negate();
+        d[1][ 9] = d[0][ 9].negate();
+        d[1][10] = d[0][10].negate();
+        d[1][11] = d[0][11].negate();
+        d[1][12] = d[0][12].negate(); // really 0
+        d[1][13] = d[0][13].negate(); // really 0
+        d[1][14] = d[0][14].negate(); // really 0
+        d[1][15] = d[0][15].negate(); // really 0
+
+        d[2][ 0] = d[0][ 0].multiply(2).subtract(1);
+        d[2][ 1] = d[0][ 1].multiply(2);
+        d[2][ 2] = d[0][ 2].multiply(2);
+        d[2][ 3] = d[0][ 3].multiply(2);
+        d[2][ 4] = d[0][ 4].multiply(2);
+        d[2][ 5] = d[0][ 5].multiply(2);
+        d[2][ 6] = d[0][ 6].multiply(2);
+        d[2][ 7] = d[0][ 7].multiply(2);
+        d[2][ 8] = d[0][ 8].multiply(2);
+        d[2][ 9] = d[0][ 9].multiply(2);
+        d[2][10] = d[0][10].multiply(2);
+        d[2][11] = d[0][11].multiply(2);
+        d[2][12] = d[0][12].multiply(2); // really 0
+        d[2][13] = d[0][13].multiply(2); // really 0
+        d[2][14] = d[0][14].multiply(2); // really 0
+        d[2][15] = d[0][15].multiply(2); // really 0
+
+        d[3][ 0] = fraction(        -17751989329.0, 2106076560.0);
+        d[3][ 1] = integrator.getField().getZero();
+        d[3][ 2] = integrator.getField().getZero();
+        d[3][ 3] = integrator.getField().getZero();
+        d[3][ 4] = integrator.getField().getZero();
+        d[3][ 5] = fraction(          4272954039.0, 7539864640.0);
+        d[3][ 6] = fraction(       -118476319744.0, 38604839385.0);
+        d[3][ 7] = fraction(        755123450731.0, 316657731600.0);
+        d[3][ 8] = fraction( 3692384461234828125.0, 1744130441634250432.0);
+        d[3][ 9] = fraction(         -4612609375.0, 5293382976.0);
+        d[3][10] = fraction(       2091772278379.0, 933644586600.0);
+        d[3][11] = fraction(          2136624137.0, 3382989120.0);
+        d[3][12] = fraction(             -126493.0, 1421424.0);
+        d[3][13] = fraction(            98350000.0, 5419179.0);
+        d[3][14] = fraction(           -18878125.0, 2053168.0);
+        d[3][15] = fraction(         -1944542619.0, 438351368.0);
+
+        d[4][ 0] = fraction(         32941697297.0, 3159114840.0);
+        d[4][ 1] = integrator.getField().getZero();
+        d[4][ 2] = integrator.getField().getZero();
+        d[4][ 3] = integrator.getField().getZero();
+        d[4][ 4] = integrator.getField().getZero();
+        d[4][ 5] = fraction(        456696183123.0, 1884966160.0);
+        d[4][ 6] = fraction(      19132610714624.0, 115814518155.0);
+        d[4][ 7] = fraction(    -177904688592943.0, 474986597400.0);
+        d[4][ 8] = fraction(-4821139941836765625.0, 218016305204281304.0);
+        d[4][ 9] = fraction(         30702015625.0, 3970037232.0);
+        d[4][10] = fraction(     -85916079474274.0, 2800933759800.0);
+        d[4][11] = fraction(         -5919468007.0, 634310460.0);
+        d[4][12] = fraction(             2479159.0, 157936.0);
+        d[4][13] = fraction(           -18750000.0, 602131.0);
+        d[4][14] = fraction(           -19203125.0, 2053168.0);
+        d[4][15] = fraction(         15700361463.0, 438351368.0);
+
+        d[5][ 0] = fraction(         12627015655.0, 631822968.0);
+        d[5][ 1] = integrator.getField().getZero();
+        d[5][ 2] = integrator.getField().getZero();
+        d[5][ 3] = integrator.getField().getZero();
+        d[5][ 4] = integrator.getField().getZero();
+        d[5][ 5] = fraction(        -72955222965.0, 188496616.0);
+        d[5][ 6] = fraction(     -13145744952320.0, 69488710893.0);
+        d[5][ 7] = fraction(      30084216194513.0, 56998391688.0);
+        d[5][ 8] = fraction( -296858761006640625.0, 25648977082856624.0);
+        d[5][ 9] = fraction(           569140625.0, 82709109.0);
+        d[5][10] = fraction(        -18684190637.0, 18672891732.0);
+        d[5][11] = fraction(            69644045.0, 89549712.0);
+        d[5][12] = fraction(           -11847025.0, 4264272.0);
+        d[5][13] = fraction(          -978650000.0, 16257537.0);
+        d[5][14] = fraction(           519371875.0, 6159504.0);
+        d[5][15] = fraction(          5256837225.0, 438351368.0);
+
+        d[6][ 0] = fraction(          -450944925.0, 17550638.0);
+        d[6][ 1] = integrator.getField().getZero();
+        d[6][ 2] = integrator.getField().getZero();
+        d[6][ 3] = integrator.getField().getZero();
+        d[6][ 4] = integrator.getField().getZero();
+        d[6][ 5] = fraction(        -14532122925.0, 94248308.0);
+        d[6][ 6] = fraction(       -595876966400.0, 2573655959.0);
+        d[6][ 7] = fraction(        188748653015.0, 527762886.0);
+        d[6][ 8] = fraction( 2545485458115234375.0, 27252038150535163.0);
+        d[6][ 9] = fraction(         -1376953125.0, 36759604.0);
+        d[6][10] = fraction(         53995596795.0, 518691437.0);
+        d[6][11] = fraction(           210311225.0, 7047894.0);
+        d[6][12] = fraction(            -1718875.0, 39484.0);
+        d[6][13] = fraction(            58000000.0, 602131.0);
+        d[6][14] = fraction(            -1546875.0, 39484.0);
+        d[6][15] = fraction(         -1262172375.0, 8429834.0);
 
     }
 
@@ -312,74 +186,6 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
 
         super(interpolator);
 
-        if (interpolator.currentState == null) {
-
-            yDotKLast = null;
-            v         = null;
-            vectorsInitialized = false;
-
-        } else {
-
-            final int dimension = interpolator.currentState.length;
-            final Field<T> field = interpolator.getGlobalPreviousState().getTime().getField();
-
-            yDotKLast = MathArrays.buildArray(field, 3, dimension);
-            for (int k = 0; k < yDotKLast.length; ++k) {
-                System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
-                                 dimension);
-            }
-
-            v = MathArrays.buildArray(field, 7, dimension);
-            for (int k = 0; k < v.length; ++k) {
-                System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
-            }
-
-            vectorsInitialized = interpolator.vectorsInitialized;
-
-        }
-
-        b_01   = interpolator.b_01;
-        b_06   = interpolator.b_06;
-        b_07   = interpolator.b_07;
-        b_08   = interpolator.b_08;
-        b_09   = interpolator.b_09;
-        b_10   = interpolator.b_10;
-        b_11   = interpolator.b_11;
-        b_12   = interpolator.b_12;
-        c14    = interpolator.c14;
-        k14_01 = interpolator.k14_01;
-        k14_06 = interpolator.k14_06;
-        k14_07 = interpolator.k14_07;
-        k14_08 = interpolator.k14_08;
-        k14_09 = interpolator.k14_09;
-        k14_10 = interpolator.k14_10;
-        k14_11 = interpolator.k14_11;
-        k14_12 = interpolator.k14_12;
-        k14_13 = interpolator.k14_13;
-        c15    = interpolator.c15;
-        k15_01 = interpolator.k15_01;
-        k15_06 = interpolator.k15_06;
-        k15_07 = interpolator.k15_07;
-        k15_08 = interpolator.k15_08;
-        k15_09 = interpolator.k15_09;
-        k15_10 = interpolator.k15_10;
-        k15_11 = interpolator.k15_11;
-        k15_12 = interpolator.k15_12;
-        k15_13 = interpolator.k15_13;
-        k15_14 = interpolator.k15_14;
-        c16    = interpolator.c16;
-        k16_01 = interpolator.k16_01;
-        k16_06 = interpolator.k16_06;
-        k16_07 = interpolator.k16_07;
-        k16_08 = interpolator.k16_08;
-        k16_09 = interpolator.k16_09;
-        k16_10 = interpolator.k16_10;
-        k16_11 = interpolator.k16_11;
-        k16_12 = interpolator.k16_12;
-        k16_13 = interpolator.k16_13;
-        k16_14 = interpolator.k16_14;
-        k16_15 = interpolator.k16_15;
-
         d = MathArrays.buildArray(integrator.getField(), 4, -1);
         for (int i = 0; i < d.length; ++i) {
             d[i] = interpolator.d[i].clone();
@@ -403,72 +209,13 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
     }
 
     /** {@inheritDoc} */
-    @Override
-    public void storeState(final FieldODEStateAndDerivative<T> state) {
-        super.storeState(state);
-        vectorsInitialized = false;
-    }
-
-    /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
     @Override
     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
                                                                                    final T time, final T theta,
                                                                                    final T oneMinusThetaH)
         throws MaxCountExceededException {
 
-        if (! vectorsInitialized) {
-
-            if (v == null) {
-                v = MathArrays.buildArray(time.getField(), 7, previousState.length);
-            }
-
-            // perform the last evaluations if they have not been done yet
-            finalizeStep();
-
-            // compute the interpolation vectors for this time step
-            for (int i = 0; i < previousState.length; ++i) {
-                final T yDot01 = yDotK[0][i];
-                final T yDot06 = yDotK[5][i];
-                final T yDot07 = yDotK[6][i];
-                final T yDot08 = yDotK[7][i];
-                final T yDot09 = yDotK[8][i];
-                final T yDot10 = yDotK[9][i];
-                final T yDot11 = yDotK[10][i];
-                final T yDot12 = yDotK[11][i];
-                final T yDot13 = yDotK[12][i];
-                final T yDot14 = yDotKLast[0][i];
-                final T yDot15 = yDotKLast[1][i];
-                final T yDot16 = yDotKLast[2][i];
-                v[0][i] =     yDot01.multiply(b_01).
-                          add(yDot06.multiply(b_06)).
-                          add(yDot07.multiply(b_07)).
-                          add(yDot08.multiply(b_08)).
-                          add(yDot09.multiply(b_09)).
-                          add(yDot10.multiply(b_10)).
-                          add(yDot11.multiply(b_11)).
-                          add(yDot12.multiply(b_12));
-                v[1][i] = yDot01.subtract(v[0][i]);
-                v[2][i] = v[0][i].subtract(v[1][i]).subtract(yDotK[12][i]);
-                for (int k = 0; k < d.length; ++k) {
-                    v[k+3][i] =     yDot01.multiply(d[k][ 0]).
-                                add(yDot06.multiply(d[k][ 1])).
-                                add(yDot07.multiply(d[k][ 2])).
-                                add(yDot08.multiply(d[k][ 3])).
-                                add(yDot09.multiply(d[k][ 4])).
-                                add(yDot10.multiply(d[k][ 5])).
-                                add(yDot11.multiply(d[k][ 6])).
-                                add(yDot12.multiply(d[k][ 7])).
-                                add(yDot13.multiply(d[k][ 8])).
-                                add(yDot14.multiply(d[k][ 9])).
-                                add(yDot15.multiply(d[k][10])).
-                                add(yDot16.multiply(d[k][11]));
-                }
-            }
-
-            vectorsInitialized = true;
-
-        }
-
         final T one      = theta.getField().getOne();
         final T eta      = one.subtract(theta);
         final T twoTheta = theta.multiply(2);
@@ -479,111 +226,32 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
         final T dot4     = theta2.multiply(theta.multiply(theta.multiply(5).subtract(8)).add(3));
         final T dot5     = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(-6).add(15)).subtract(12)).add(3));
         final T dot6     = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-7).add(18)).subtract(15)).add(4)));
-        final T[] interpolatedState       = MathArrays.buildArray(theta.getField(), previousState.length);
-        final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
-
-        if ((previousState != null) && (theta.getReal() <= 0.5)) {
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                interpolatedState[i] = previousState[i].
-                                       add(theta.multiply(h.multiply(v[0][i].
-                                                          add(eta.multiply(v[1][i].
-                                                                           add(theta.multiply(v[2][i].
-                                                                                              add(eta.multiply(v[3][i].
-                                                                                                               add(theta.multiply(v[4][i].
-                                                                                                                   add(eta.multiply(v[5][i].
-                                                                                                                                    add(theta.multiply(v[6][i])))))))))))))));
-                interpolatedDerivatives[i] =                   v[0][i].
-                                             add(dot1.multiply(v[1][i])).
-                                             add(dot2.multiply(v[2][i])).
-                                             add(dot3.multiply(v[3][i])).
-                                             add(dot4.multiply(v[4][i])).
-                                             add(dot5.multiply(v[5][i])).
-                                             add(dot6.multiply(v[6][i]));
-            }
+        final T[] interpolatedState;
+        final T[] interpolatedDerivatives;
+
+        if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+            final T f0 = theta.multiply(h);
+            final T f1 = f0.multiply(eta);
+            final T f2 = f1.multiply(theta);
+            final T f3 = f2.multiply(eta);
+            final T f4 = f3.multiply(theta);
+            final T f5 = f4.multiply(eta);
+            final T f6 = f5.multiply(theta);
+            interpolatedState       = previousStateLinearCombination(f0, f1, f2, f3, f4, f5, f6);
+            interpolatedDerivatives = derivativeLinearCombination(one, dot1, dot2, dot3, dot4, dot5, dot6);
         } else {
-            for (int i = 0; i < interpolatedState.length; ++i) {
-                interpolatedState[i] = currentState[i].
-                                       subtract(oneMinusThetaH.multiply(v[0][i].
-                                                                        subtract(theta.multiply(v[1][i].
-                                                                                                add(theta.multiply(v[2][i].
-                                                                                                                   add(eta.multiply(v[3][i].
-                                                                                                                                    add(theta.multiply(v[4][i].
-                                                                                                                                                       add(eta.multiply(v[5][i].
-                                                                                                                                                                        add(theta.multiply(v[6][i]))))))))))))));
-                interpolatedDerivatives[i] =                   v[0][i].
-                                             add(dot1.multiply(v[1][i])).
-                                             add(dot2.multiply(v[2][i])).
-                                             add(dot3.multiply(v[3][i])).
-                                             add(dot4.multiply(v[4][i])).
-                                             add(dot5.multiply(v[5][i])).
-                                             add(dot6.multiply(v[6][i]));
-            }
+            final T f0 = oneMinusThetaH;
+            final T f1 = f0.multiply(theta).negate();
+            final T f2 = f1.multiply(theta);
+            final T f3 = f2.multiply(eta);
+            final T f4 = f3.multiply(theta);
+            final T f5 = f4.multiply(eta);
+            final T f6 = f5.multiply(theta);
+            interpolatedState       = currentStateLinearCombination(f0, f1, f2, f3, f4, f5, f6);
+            interpolatedDerivatives = derivativeLinearCombination(one, dot1, dot2, dot3, dot4, dot5, dot6);
         }
 
-        return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
-
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    protected void doFinalize() throws MaxCountExceededException {
-
-        if (currentState == null) {
-            // we are finalizing an uninitialized instance
-            return;
-        }
-
-        T s;
-        final T   pT   = getGlobalPreviousState().getTime();
-        final T[] yTmp = MathArrays.buildArray(pT.getField(), currentState.length);
-
-        // k14
-        for (int j = 0; j < currentState.length; ++j) {
-            s =     yDotK[ 0][j].multiply(k14_01).
-                add(yDotK[ 5][j].multiply(k14_06)).
-                add(yDotK[ 6][j].multiply(k14_07)).
-                add(yDotK[ 7][j].multiply(k14_08)).
-                add(yDotK[ 8][j].multiply(k14_09)).
-                add(yDotK[ 9][j].multiply(k14_10)).
-                add(yDotK[10][j].multiply(k14_11)).
-                add(yDotK[11][j].multiply(k14_12)).
-                add(yDotK[12][j].multiply(k14_13));
-            yTmp[j] = currentState[j].add(h.multiply(s));
-        }
-        yDotKLast[0] = integrator.computeDerivatives(pT.add(h.multiply(c14)), yTmp);
-
-        // k15
-        for (int j = 0; j < currentState.length; ++j) {
-            s =     yDotK[ 0][j].multiply(k15_01).
-                add(yDotK[ 5][j].multiply(k15_06)).
-                add(yDotK[ 6][j].multiply(k15_07)).
-                add(yDotK[ 7][j].multiply(k15_08)).
-                add(yDotK[ 8][j].multiply(k15_09)).
-                add(yDotK[ 9][j].multiply(k15_10)).
-                add(yDotK[10][j].multiply(k15_11)).
-                add(yDotK[11][j].multiply(k15_12)).
-                add(yDotK[12][j].multiply(k15_13)).
-                add(yDotKLast[0][j].multiply(k15_14));
-            yTmp[j] = currentState[j].add(h.multiply(s));
-        }
-        yDotKLast[1] = integrator.computeDerivatives(pT.add(h.multiply(c15)), yTmp);
-
-        // k16
-        for (int j = 0; j < currentState.length; ++j) {
-            s =     yDotK[ 0][j].multiply(k16_01).
-                add(yDotK[ 5][j].multiply(k16_06)).
-                add(yDotK[ 6][j].multiply(k16_07)).
-                add(yDotK[ 7][j].multiply(k16_08)).
-                add(yDotK[ 8][j].multiply(k16_09)).
-                add(yDotK[ 9][j].multiply(k16_10)).
-                add(yDotK[10][j].multiply(k16_11)).
-                add(yDotK[11][j].multiply(k16_12)).
-                add(yDotK[12][j].multiply(k16_13)).
-                add(yDotKLast[0][j].multiply(k16_14)).
-                add(yDotKLast[0][j].multiply(k16_15));
-            yTmp[j] = currentState[j].add(h.multiply(s));
-        }
-        yDotKLast[2] = integrator.computeDerivatives(pT.add(h.multiply(c16)), yTmp);
+        return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
index 8da12a9..5f648f8 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
@@ -23,7 +23,6 @@ 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.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.ode.FieldExpandableODE;
 import org.apache.commons.math3.ode.FieldODEState;
@@ -71,8 +70,8 @@ import org.apache.commons.math3.util.MathUtils;
 public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
     extends AdaptiveStepsizeFieldIntegrator<T> {
 
-    /** Indicator for <i>fsal</i> methods. */
-    private final boolean fsal;
+    /** Index of the pre-computed derivative for <i>fsal</i> methods. */
+    private final int fsal;
 
     /** Time steps from Butcher array (without the first zero). */
     private final T[] c;
@@ -98,7 +97,8 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
     /** Build a Runge-Kutta integrator with the given Butcher array.
      * @param field field to which the time and state vector elements belong
      * @param name name of the method
-     * @param fsal indicate that the method is an <i>fsal</i>
+     * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
+     * or -1 if method is not <i>fsal</i>
      * @param minStep minimal step (sign is irrelevant, regardless of
      * integration direction, forward or backward), the last step can
      * be smaller than this
@@ -108,7 +108,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
      * @param scalAbsoluteTolerance allowed absolute error
      * @param scalRelativeTolerance allowed relative error
      */
-    protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final boolean fsal,
+    protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
                                                 final double minStep, final double maxStep,
                                                 final double scalAbsoluteTolerance,
                                                 final double scalRelativeTolerance) {
@@ -132,7 +132,8 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
     /** Build a Runge-Kutta integrator with the given Butcher array.
      * @param field field to which the time and state vector elements belong
      * @param name name of the method
-     * @param fsal indicate that the method is an <i>fsal</i>
+     * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
+     * or -1 if method is not <i>fsal</i>
      * @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
@@ -140,7 +141,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
      * @param vecAbsoluteTolerance allowed absolute error
      * @param vecRelativeTolerance allowed relative error
      */
-    protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final boolean fsal,
+    protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
                                                 final double   minStep, final double maxStep,
                                                 final double[] vecAbsoluteTolerance,
                                                 final double[] vecRelativeTolerance) {
@@ -195,18 +196,11 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
     protected abstract T[] getB();
 
     /** Create an interpolator.
-     * @param rkIntegrator integrator being used
-     * @param y reference to the integrator array holding the state at
-     * the end of the step
-     * @param yDotArray reference to the integrator array holding all the
-     * intermediate slopes
      * @param forward integration direction indicator
      * @param mapper equations mapper for the all equations
      * @return external weights for the high order method from Butcher array
      */
-    protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(AbstractFieldIntegrator<T> rkIntegrator,
-                                                                             T[] y, T[][] yDotArray,
-                                                                             boolean forward,
+    protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward,
                                                                              FieldEquationsMapper<T> mapper);
     /** Get the order of the method.
      * @return order of the method
@@ -248,7 +242,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
 
         // set up an interpolator sharing the integrator arrays
         final RungeKuttaFieldStepInterpolator<T> interpolator =
-                        createInterpolator(this, y0, yDotK, forward, equations.getMapper());
+                        createInterpolator(forward, equations.getMapper());
         interpolator.storeState(stepStart);
 
         // set up integration control objects
@@ -328,11 +322,12 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
                 }
 
             }
-            final T stepEnd   = stepStart.getTime().add(stepSize);
-            final T[] yDotTmp = fsal ? yDotK[stages - 1] : computeDerivatives(stepEnd, yTmp);
+            final T   stepEnd = stepStart.getTime().add(stepSize);
+            final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
             final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
 
             // local error is small enough: accept the step, trigger events and step handlers
+            interpolator.setSlopes(yDotK);
             interpolator.storeState(stateTmp);
             System.arraycopy(yTmp, 0, y, 0, y0.length);
             stepStart = acceptStep(interpolator, finalTime);

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
index 7fb32a1..4c58e09 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
 
 import org.apache.commons.math3.Field;
 import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.util.MathArrays;
 
@@ -86,10 +85,8 @@ public class EulerFieldIntegrator<T extends RealFieldElement<T>> extends RungeKu
     /** {@inheritDoc} */
     @Override
     protected EulerFieldStepInterpolator<T>
-        createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
-                           final T[][] yDotArray, final boolean forward,
-                           final FieldEquationsMapper<T> mapper) {
-        return new EulerFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+        createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+        return new EulerFieldStepInterpolator<T>(this, forward, mapper);
     }
 
 }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
index dba2312..9bff5f2 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
 import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
 
 /**
  * This class implements a linear interpolator for step.
@@ -52,17 +51,13 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
 
     /** Simple constructor.
      * @param rkIntegrator integrator being used
-     * @param y reference to the integrator array holding the state at
-     * the end of the step
-     * @param yDotArray reference to the integrator array holding all the
-     * intermediate slopes
      * @param forward integration direction indicator
      * @param mapper equations mapper for the all equations
      */
     EulerFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
-                               final T[] y, final T[][] yDotArray, final boolean forward,
+                               final boolean forward,
                                final FieldEquationsMapper<T> mapper) {
-        super(rkIntegrator, y, yDotArray, forward, mapper);
+        super(rkIntegrator, forward, mapper);
     }
 
     /** Copy constructor.
@@ -81,22 +76,23 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
     }
 
     /** {@inheritDoc} */
+    @SuppressWarnings("unchecked")
     @Override
     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
                                                                                    final T time, final T theta,
                                                                                    final T oneMinusThetaH) {
-        final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
-        if ((previousState != null) && (theta.getReal() <= 0.5)) {
-            for (int i = 0; i < previousState.length; ++i) {
-                interpolatedState[i] = previousState[i].add(theta.multiply(h).multiply(yDotK[0][i]));
-            }
+        final T[] interpolatedState;
+        final T[] interpolatedDerivatives;
+        if ((getGlobalPreviousState() != null) && (theta.getReal() <= 0.5)) {
+            interpolatedState       = previousStateLinearCombination(theta.multiply(h));
+            interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
         } else {
-            for (int i = 0; i < previousState.length; ++i) {
-                interpolatedState[i] = currentState[i].subtract(oneMinusThetaH.multiply(yDotK[0][i]));
-            }
+            interpolatedState       = currentStateLinearCombination(oneMinusThetaH.negate());
+            interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
         }
 
-        return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+        return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
+
     }
 
 }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
index d8460c2..659eee8 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
 
 import org.apache.commons.math3.Field;
 import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
 import org.apache.commons.math3.ode.FieldEquationsMapper;
 import org.apache.commons.math3.util.MathArrays;
 
@@ -111,10 +110,8 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
     /** {@inheritDoc} */
     @Override
     protected GillFieldStepInterpolator<T>
-        createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
-                           final T[][] yDotArray, final boolean forward,
-                           final FieldEquationsMapper<T> mapper) {
-        return new GillFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+        createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+        return new GillFieldStepInterpolator<T>(this, forward, mapper);
     }
 
 }


Mime
View raw message