http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java index 8a6b38e..8990dc5 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java @@ -22,7 +22,6 @@ 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.FastMath; -import org.apache.commons.math3.util.MathArrays; /** * This class implements a step interpolator for the Gill fourth @@ -68,17 +67,13 @@ class GillFieldStepInterpolator<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 */ GillFieldStepInterpolator(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. @@ -98,6 +93,7 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>> /** {@inheritDoc} */ + @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, @@ -111,48 +107,30 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>> final T coeffDot2 = cDot23.multiply(ONE_MINUS_INV_SQRT_2); final T coeffDot3 = cDot23.multiply(ONE_PLUS_INV_SQRT_2); final T coeffDot4 = theta.multiply(twoTheta.subtract(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)) { - final T s = theta.multiply(h).divide(6.0); - final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2)); - final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6)); - final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2); - final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2); - 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 yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - interpolatedState[i] = previousState[i]. - add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)). - add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4)); - } + final T[] interpolatedState; + final T[] interpolatedDerivatives; + + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { + final T s = theta.multiply(h).divide(6.0); + final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2)); + final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6)); + final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2); + final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2); + final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3))); + interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4); } else { - final T s = oneMinusThetaH.divide(6.0); - final T c23 = s .multiply(twoTheta.add(2).subtract(fourTheta2)); + final T s = oneMinusThetaH.divide(-6.0); + final T c23 = s.multiply(twoTheta.add(2).subtract(fourTheta2)); final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(5)).add(1)); final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2); final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2); final T coeff4 = s.multiply(fourTheta2.add(theta).add(1)); - for (int i = 0; i < interpolatedState.length; ++i) { - final T yDot1 = yDotK[0][i]; - final T yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - interpolatedState[i] = currentState[i]. - subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)). - subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4)); - } + interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, 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/HighamHall54FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java index 092393a..8ee90b7 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.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; @@ -64,7 +63,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { - super(field, METHOD_NAME, false, + super(field, METHOD_NAME, -1, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); e = MathArrays.buildArray(field, 7); e[0] = fraction(-1, 20); @@ -92,7 +91,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { - super(field, METHOD_NAME, false, + super(field, METHOD_NAME, -1, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); e = MathArrays.buildArray(field, 7); e[0] = fraction(-1, 20); @@ -165,10 +164,8 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected HighamHall54FieldStepInterpolator<T> - createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y, - final T[][] yDotArray, final boolean forward, - final FieldEquationsMapper<T> mapper) { - return new HighamHall54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper); + createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { + return new HighamHall54FieldStepInterpolator<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/HighamHall54FieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java index 3ee52a9..3f8499c 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.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 @@ -38,17 +37,13 @@ class HighamHall54FieldStepInterpolator<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 */ HighamHall54FieldStepInterpolator(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. @@ -68,72 +63,44 @@ class HighamHall54FieldStepInterpolator<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 bDot0 = theta.multiply(theta.multiply(theta.multiply( -10.0 ).add( 16.0 )).add(-15.0 / 2.0)).add(1); + final T bDot1 = time.getField().getZero(); final T bDot2 = theta.multiply(theta.multiply(theta.multiply( 135.0 / 2.0).add(-729.0 / 8.0)).add(459.0 / 16.0)); final T bDot3 = theta.multiply(theta.multiply(theta.multiply(-120.0 ).add( 152.0 )).add(-44.0 )); final T bDot4 = theta.multiply(theta.multiply(theta.multiply( 125.0 / 2.0).add(-625.0 / 8.0)).add(375.0 / 16.0)); final T bDot5 = theta.multiply( 5.0 / 8.0).multiply(theta.multiply(2).subtract(1)); - 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)) { + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { final T hTheta = h.multiply(theta); final T b0 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply( -5.0 / 2.0).add( 16.0 / 3.0)).add(-15.0 / 4.0)).add(1)); + final T b1 = time.getField().getZero(); final T b2 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0))); final T b3 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 ))); final T b4 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0))); final T b5 = hTheta.multiply(theta.multiply(theta.multiply( 5.0 / 12.0)).add( -5.0 / 16.0)); - for (int i = 0; i < interpolatedState.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]; - interpolatedState[i] = previousState[i]. - add(b0.multiply(yDot0)). - add(b2.multiply(yDot2)). - add(b3.multiply(yDot3)). - add(b4.multiply(yDot4)). - add(b5.multiply(yDot5)); - interpolatedDerivatives[i] = bDot0.multiply(yDot0). - add(bDot2.multiply(yDot2)). - add(bDot3.multiply(yDot3)). - add(bDot4.multiply(yDot4)). - add(bDot5.multiply(yDot5)); - } + interpolatedState = previousStateLinearCombination(b0, b1, b2, b3, b4, b5); + interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5); } else { final T theta2 = theta.multiply(theta); final T b0 = h.multiply( theta.multiply(theta.multiply(theta.multiply(theta.multiply(-5.0 / 2.0).add( 16.0 / 3.0)).add( -15.0 / 4.0)).add( 1.0 )).add( -1.0 / 12.0)); + final T b1 = time.getField().getZero(); final T b2 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 135.0 / 8.0 ).add(-243.0 / 8.0)).add(459.0 / 32.0)).add( -27.0 / 32.0)); final T b3 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( -30.0 ).add( 152.0 / 3.0)).add(-22.0 )).add( 4.0 / 3.0)); final T b4 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 125.0 / 8.0 ).add(-625.0 / 24.0)).add(375.0 / 32.0)).add(-125.0 / 96.0)); final T b5 = h.multiply(theta2.multiply(theta.multiply( 5.0 / 12.0 ).add(-5.0 / 16.0)).add( -5.0 / 48.0)); - for (int i = 0; i < interpolatedState.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]; - interpolatedState[i] = currentState[i]. - add(b0.multiply(yDot0)). - add(b2.multiply(yDot2)). - add(b3.multiply(yDot3)). - add(b4.multiply(yDot4)). - add(b5.multiply(yDot5)); - interpolatedDerivatives[i] = bDot0.multiply(yDot0). - add(bDot2.multiply(yDot2)). - add(bDot3.multiply(yDot3)). - add(bDot4.multiply(yDot4)). - add(bDot5.multiply(yDot5)); - } + interpolatedState = currentStateLinearCombination(b0, b1, b2, b3, b4, b5); + interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5); } - 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/LutherFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java index 363e042..c36e5dd 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.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; @@ -138,10 +137,8 @@ public class LutherFieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected LutherFieldStepInterpolator<T> - createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y, - final T[][] yDotArray, final boolean forward, - final FieldEquationsMapper<T> mapper) { - return new LutherFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper); + createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { + return new LutherFieldStepInterpolator<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/LutherFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java index 8a0d407..6f643d5 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.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 @@ -83,17 +82,13 @@ class LutherFieldStepInterpolator<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 */ LutherFieldStepInterpolator(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); final T q = rkIntegrator.getField().getOne().multiply(21).sqrt(); c5a = q.multiply( -49).add( -49); c5b = q.multiply( 287).add( 392); @@ -143,6 +138,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> /** {@inheritDoc} */ + @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, @@ -192,86 +188,42 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> // At the end, we get the b_i as polynomials in theta. final T coeffDot1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 ).add( -47 )).add( 36 )).add( -54 / 5.0)).add(1); - // not really needed as it is zero: final T coeffDot2 = theta.getField().getZero(); + final T coeffDot2 = theta.getField().getZero(); final T coeffDot3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 ).add(-608 / 3.0)).add( 320 / 3.0 )).add(-208 / 15.0)); final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0)); final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150))); final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150))); final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3 )).add( -3 )).add( 3 / 5.0); - 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)) { - - final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1); - // not really needed as it is zero: final T coeff2 = theta.getField().getZero(); - final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0)); - final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0)); - final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))); - final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))); - final T coeff7 = theta.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0); - for (int i = 0; i < interpolatedState.length; ++i) { - final T yDot1 = yDotK[0][i]; - // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - final T yDot5 = yDotK[4][i]; - final T yDot6 = yDotK[5][i]; - final T yDot7 = yDotK[6][i]; - interpolatedState[i] = previousState[i]. - add(theta.multiply(h). - multiply( coeff1.multiply(yDot1). - // not really needed as it is zero: add(coeff2.multiply(yDot2)). - add(coeff3.multiply(yDot3)). - add(coeff4.multiply(yDot4)). - add(coeff5.multiply(yDot5)). - add(coeff6.multiply(yDot6)). - add(coeff7.multiply(yDot7)))); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1). - // not really needed as it is zero: add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)). - add(coeffDot4.multiply(yDot4)). - add(coeffDot5.multiply(yDot5)). - add(coeffDot6.multiply(yDot6)). - add(coeffDot7.multiply(yDot7)); - } + final T[] interpolatedState; + final T[] interpolatedDerivatives; + + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { + + final T s = theta.multiply(theta.multiply(h)); + final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1); + final T coeff2 = s.getField().getZero(); + final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0)); + final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0)); + final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))); + final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))); + final T coeff7 = s.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0); + interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7); } else { - final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0); - // not really needed as it is zero: final T coeff2 = theta.getField().getZero(); - final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0); - final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0))); - final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); - final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); - final T coeff7 = theta.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0); - for (int i = 0; i < interpolatedState.length; ++i) { - final T yDot1 = yDotK[0][i]; - // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - final T yDot5 = yDotK[4][i]; - final T yDot6 = yDotK[5][i]; - final T yDot7 = yDotK[6][i]; - interpolatedState[i] = currentState[i]. - add(oneMinusThetaH. - multiply( coeff1.multiply(yDot1). - // not really needed as it is zero: add(coeff2.multiply(yDot2)). - add(coeff3.multiply(yDot3)). - add(coeff4.multiply(yDot4)). - add(coeff5.multiply(yDot5)). - add(coeff6.multiply(yDot6)). - add(coeff7.multiply(yDot7)))); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1). - // not really needed as it is zero: add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)). - add(coeffDot4.multiply(yDot4)). - add(coeffDot5.multiply(yDot5)). - add(coeffDot6.multiply(yDot6)). - add(coeffDot7.multiply(yDot7)); - } + final T s = oneMinusThetaH.multiply(theta); + final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0); + final T coeff2 = s.getField().getZero(); + final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0); + final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0))); + final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); + final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); + final T coeff7 = s.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0); + interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7); } - 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/MidpointFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java index f74fdba..3cc4c71 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.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 MidpointFieldIntegrator<T extends RealFieldElement<T>> extends Rung /** {@inheritDoc} */ @Override protected MidpointFieldStepInterpolator<T> - createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y, - final T[][] yDotArray, final boolean forward, - final FieldEquationsMapper<T> mapper) { - return new MidpointFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper); + createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { + return new MidpointFieldStepInterpolator<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/MidpointFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java index 6569c4c..943533a 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.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 second order @@ -54,17 +53,13 @@ class MidpointFieldStepInterpolator<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 */ MidpointFieldStepInterpolator(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. @@ -84,37 +79,30 @@ class MidpointFieldStepInterpolator<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 coeffDot2 = theta.multiply(2); - final T coeffDot1 = time.getField().getOne().subtract(coeffDot2); - final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length); - final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length); + final T coeffDot2 = theta.multiply(2); + final T coeffDot1 = time.getField().getOne().subtract(coeffDot2); + final T[] interpolatedState; + final T[] interpolatedDerivatives; - if ((previousState != null) && (theta.getReal() <= 0.5)) { + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { final T coeff1 = theta.multiply(oneMinusThetaH); final T coeff2 = theta.multiply(theta).multiply(h); - for (int i = 0; i < previousState.length; ++i) { - final T yDot1 = yDotK[0][i]; - final T yDot2 = yDotK[1][i]; - interpolatedState[i] = previousState[i].add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)); - } + interpolatedState = previousStateLinearCombination(coeff1, coeff2); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2); } else { final T coeff1 = oneMinusThetaH.multiply(theta); - final T coeff2 = oneMinusThetaH.multiply(theta.add(1)); - for (int i = 0; i < previousState.length; ++i) { - final T yDot1 = yDotK[0][i]; - final T yDot2 = yDotK[1][i]; - interpolatedState[i] = currentState[i].add(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)); - } + final T coeff2 = oneMinusThetaH.multiply(theta.add(1)).negate(); + interpolatedState = currentStateLinearCombination(coeff1, coeff2); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2); } - 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/RungeKuttaFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java index 61750f6..a511ac8 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java @@ -112,18 +112,11 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> 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); /** {@inheritDoc} */ @@ -147,7 +140,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> // 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 @@ -172,6 +165,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> interpolator.shift(); // first stage + y = equations.getMapper().mapState(stepStart); yDotK[0] = equations.getMapper().mapDerivative(stepStart); // next stages @@ -202,6 +196,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>> final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp); // discrete events handling + 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/RungeKuttaFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java index fcc0396..9a86d82 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java @@ -36,31 +36,23 @@ import org.apache.commons.math3.util.MathArrays; abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> extends AbstractFieldStepInterpolator<T> { - /** Previous state. */ - protected T[] previousState; - - /** Slopes at the intermediate points */ - protected T[][] yDotK; - /** Reference to the integrator. */ protected AbstractFieldIntegrator<T> integrator; + /** Slopes at the intermediate points. */ + private T[][] yDotK; + /** 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 */ protected RungeKuttaFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, - final T[] y, final T[][] yDotArray, final boolean forward, + final boolean forward, final FieldEquationsMapper<T> mapper) { - super(y, forward, mapper); - this.previousState = null; - this.yDotK = yDotArray; - this.integrator = rkIntegrator; + super(forward, mapper); + this.yDotK = null; + this.integrator = rkIntegrator; } /** Copy constructor. @@ -84,18 +76,14 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> super(interpolator); - if (interpolator.currentState != null) { - - previousState = interpolator.previousState.clone(); - + if (yDotK != null) { yDotK = MathArrays.buildArray(interpolator.integrator.getField(), - interpolator.yDotK.length, interpolator.yDotK[0].length); - for (int k = 0; k < interpolator.yDotK.length; ++k) { - System.arraycopy(interpolator.yDotK[k], 0, yDotK[k], 0, interpolator.yDotK[k].length); + interpolator.yDotK.length, -1); + for (int k = 0; k < yDotK.length; ++k) { + yDotK[k] = interpolator.yDotK[k].clone(); } } else { - previousState = null; yDotK = null; } @@ -105,11 +93,56 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> } - /** {@inheritDoc} */ - @Override - public void shift() { - previousState = currentState.clone(); - super.shift(); + /** Store the slopes at the intermediate points. + * @param slopes slopes at the intermediate points + */ + void setSlopes(final T[][] slopes) { + this.yDotK = slopes.clone(); + } + + /** Compute a state by linear combination added to previous state. + * @param coefficients coefficients to apply to the method staged derivatives + * @return combined state + */ + @SafeVarargs + protected final T[] previousStateLinearCombination(final T ... coefficients) { + return combine(getPreviousState().getState(), + coefficients); + } + + /** Compute a state by linear combination added to current state. + * @param coefficients coefficients to apply to the method staged derivatives + * @return combined state + */ + @SuppressWarnings("unchecked") + protected T[] currentStateLinearCombination(final T ... coefficients) { + return combine(getCurrentState().getState(), + coefficients); + } + + /** Compute a state derivative by linear combination. + * @param coefficients coefficients to apply to the method staged derivatives + * @return combined state + */ + @SuppressWarnings("unchecked") + protected T[] derivativeLinearCombination(final T ... coefficients) { + return combine(MathArrays.buildArray(integrator.getField(), yDotK[0].length), + coefficients); + } + + /** Linearly combine arrays. + * @param a array to add to + * @param coefficients coefficients to apply to the method staged derivatives + * @return a itself, as a conveniency for fluent API + */ + @SuppressWarnings("unchecked") + private T[] combine(final T[] a, final T ... coefficients) { + for (int i = 0; i < a.length; ++i) { + for (int k = 0; k < coefficients.length; ++k) { + a[i] = a[i].add(coefficients[k].multiply(yDotK[k][i])); + } + } + return a; } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java index 6ba4b72..e0d5044 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.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; @@ -100,10 +99,8 @@ public class ThreeEighthesFieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected ThreeEighthesFieldStepInterpolator<T> - createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y, - final T[][] yDotArray, final boolean forward, - final FieldEquationsMapper<T> mapper) { - return new ThreeEighthesFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper); + createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { + return new ThreeEighthesFieldStepInterpolator<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/ThreeEighthesFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java index 3de88f6..39cdc1c 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.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 3/8 fourth @@ -64,17 +63,13 @@ class ThreeEighthesFieldStepInterpolator<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 */ ThreeEighthesFieldStepInterpolator(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. @@ -94,6 +89,7 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>> /** {@inheritDoc} */ + @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, @@ -103,51 +99,31 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>> final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1); final T coeffDot2 = coeffDot3.multiply(theta.multiply(-6).add(5)); final T coeffDot4 = coeffDot3.multiply(theta.multiply(2).subtract(1)); - 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)) { + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { final T s = theta.multiply(h).divide(8); final T fourTheta2 = theta.multiply(theta).multiply(4); final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(15)).add(8)); final T coeff2 = s.multiply(theta.multiply(5).subtract(fourTheta2)).multiply(3); final T coeff3 = s.multiply(theta).multiply(3); 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 yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - interpolatedState[i] = previousState[i]. - add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)). - add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4)); - - } + interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4); } else { - final T s = oneMinusThetaH.divide(8); + final T s = oneMinusThetaH.divide(-8); final T fourTheta2 = theta.multiply(theta).multiply(4); final T thetaPlus1 = theta.add(1); final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(7)).add(1)); final T coeff2 = s.multiply(thetaPlus1.subtract(fourTheta2)).multiply(3); final T coeff3 = s.multiply(thetaPlus1).multiply(3); final T coeff4 = s.multiply(thetaPlus1.add(fourTheta2)); - for (int i = 0; i < interpolatedState.length; ++i) { - final T yDot1 = yDotK[0][i]; - final T yDot2 = yDotK[1][i]; - final T yDot3 = yDotK[2][i]; - final T yDot4 = yDotK[3][i]; - interpolatedState[i] = currentState[i]. - subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)). - subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4)); - interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)). - add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4)); - - } + interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4); + interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, 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/sampling/AbstractFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java index b450a77..e9d6f29 100644 --- a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java @@ -43,9 +43,6 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T /** Current time step. */ protected T h; - /** Current state. */ - protected T[] currentState; - /** Global previous state. */ private FieldODEStateAndDerivative<T> globalPreviousState; @@ -68,18 +65,16 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T private FieldEquationsMapper<T> mapper; /** Simple constructor. - * @param y reference to the integrator array holding the state at the end of the step * @param isForward integration direction indicator * @param equationsMapper mapper for ODE equations primary and secondary components */ - protected AbstractFieldStepInterpolator(final T[] y, final boolean isForward, + protected AbstractFieldStepInterpolator(final boolean isForward, final FieldEquationsMapper<T> equationsMapper) { globalPreviousState = null; globalCurrentState = null; softPreviousState = null; softCurrentState = null; h = null; - currentState = y; finalized = false; this.forward = isForward; this.mapper = equationsMapper; @@ -109,17 +104,9 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T softPreviousState = interpolator.softPreviousState; softCurrentState = interpolator.softCurrentState; h = interpolator.h; - - if (interpolator.currentState == null) { - currentState = null; - mapper = null; - } else { - currentState = interpolator.currentState.clone(); - } - - finalized = interpolator.finalized; - forward = interpolator.forward; - mapper = interpolator.mapper; + finalized = interpolator.finalized; + forward = interpolator.forward; + mapper = interpolator.mapper; } http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java index f79fad0..278276b 100644 --- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java @@ -66,9 +66,10 @@ public class TestFieldProblem1<T extends RealFieldElement<T>> @Override public T[] computeTheoreticalState(T t) { - final T[] y0 = getInitialState(); + final FieldODEState<T> s0 = getInitialState(); + final T[] y0 = s0.getState(); final T[] y = MathArrays.buildArray(getField(), getDimension()); - T c = getInitialTime().subtract(t).exp(); + T c = s0.getTime().subtract(t).exp(); for (int i = 0; i < getDimension(); ++i) { y[i] = c.multiply(y0[i]); } http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java index f84781a..80f0ea2 100644 --- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java @@ -35,7 +35,7 @@ public class TestFieldProblem5<T extends RealFieldElement<T>> */ public TestFieldProblem5(Field<T> field) { super(field); - setFinalConditions(getInitialTime().multiply(2).subtract(getFinalTime())); + setFinalConditions(getInitialState().getTime().multiply(2).subtract(getFinalTime())); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java index 09a10d8..02ff65f 100644 --- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java @@ -109,20 +109,12 @@ public abstract class TestFieldProblemAbstract<T extends RealFieldElement<T>> return n; } - /** - * Get the initial time. - * @return initial time - */ - public T getInitialTime() { - return t0; - } - - /** - * Get the initial state vector. - * @return initial state vector + /** + * Get the initial state. + * @return initial state */ - public T[] getInitialState() { - return y0; + public FieldODEState<T> getInitialState() { + return new FieldODEState<T>(t0, y0); } /** http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java index d4e20e6..101f4fa 100644 --- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java +++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java @@ -74,7 +74,7 @@ public class TestFieldProblemHandler<T extends RealFieldElement<T>> public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) throws MaxCountExceededException { T start = integrator.getCurrentStepStart().getTime(); - if (start.subtract(problem.getInitialTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) { + if (start.subtract(problem.getInitialState().getTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) { // multistep integrators do not handle the first steps themselves // so we have to make sure the integrator we look at has really started its work if (expectedStepStart != null) { http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java index 54f0330..9ec73d2 100644 --- a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java +++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java @@ -58,8 +58,8 @@ public class EulerIntegratorTest { MaxCountExceededException, NoBracketingException { for (TestProblemAbstract pb : new TestProblemAbstract[] { - new TestProblem1(), new TestProblem2(), new TestProblem3(), - new TestProblem4(), new TestProblem5(), new TestProblem6() + //new TestProblem1(), new TestProblem2(), new TestProblem3(), + new TestProblem4()//, new TestProblem5(), new TestProblem6() }) { double previousValueError = Double.NaN; @@ -94,6 +94,7 @@ public class EulerIntegratorTest { } previousTimeError = timeError; + System.exit(1); } }
