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/MATH_3_X Commit: 07e2b9447b5950721f39b7d027a0ca95d0cd950e Parents: 8180c9e Author: Luc Maisonobe <l...@apache.org> Authored: Tue Dec 1 12:22:33 2015 +0100 Committer: Luc Maisonobe <l...@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); } }
