Use immutable step interpolators.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/771eb6a6 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/771eb6a6 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/771eb6a6 Branch: refs/heads/master Commit: 771eb6a606491375a61fee1aa537025edafacc34 Parents: 0ddec29 Author: Luc Maisonobe <[email protected]> Authored: Wed Jan 6 13:22:39 2016 +0100 Committer: Luc Maisonobe <[email protected]> Committed: Wed Jan 6 13:22:39 2016 +0100 ---------------------------------------------------------------------- .../math4/ode/AbstractFieldIntegrator.java | 15 +- .../math4/ode/ContinuousOutputFieldModel.java | 6 +- .../ClassicalRungeKuttaFieldIntegrator.java | 11 +- ...lassicalRungeKuttaFieldStepInterpolator.java | 42 ++-- .../DormandPrince54FieldIntegrator.java | 10 +- .../DormandPrince54FieldStepInterpolator.java | 64 +++--- .../DormandPrince853FieldIntegrator.java | 10 +- .../DormandPrince853FieldStepInterpolator.java | 227 ++++++++++--------- .../EmbeddedRungeKuttaFieldIntegrator.java | 22 +- .../ode/nonstiff/EulerFieldIntegrator.java | 11 +- .../nonstiff/EulerFieldStepInterpolator.java | 46 ++-- .../math4/ode/nonstiff/GillFieldIntegrator.java | 11 +- .../ode/nonstiff/GillFieldStepInterpolator.java | 45 ++-- .../nonstiff/HighamHall54FieldIntegrator.java | 10 +- .../HighamHall54FieldStepInterpolator.java | 59 +++-- .../ode/nonstiff/LutherFieldIntegrator.java | 11 +- .../nonstiff/LutherFieldStepInterpolator.java | 62 +++-- .../ode/nonstiff/MidpointFieldIntegrator.java | 11 +- .../nonstiff/MidpointFieldStepInterpolator.java | 47 ++-- .../ode/nonstiff/RungeKuttaFieldIntegrator.java | 22 +- .../RungeKuttaFieldStepInterpolator.java | 79 ++++--- .../nonstiff/ThreeEighthesFieldIntegrator.java | 11 +- .../ThreeEighthesFieldStepInterpolator.java | 41 ++-- .../sampling/AbstractFieldStepInterpolator.java | 145 +++++------- .../ode/sampling/FieldStepInterpolator.java | 11 - .../ode/ContinuousOutputFieldModelTest.java | 219 ++++++++++++++++++ ...ractRungeKuttaFieldStepInterpolatorTest.java | 71 ++---- ...sicalRungKuttaFieldStepInterpolatorTest.java | 20 +- ...ormandPrince54FieldStepInterpolatorTest.java | 20 +- ...rmandPrince853FieldStepInterpolatorTest.java | 18 +- .../EulerFieldStepInterpolatorTest.java | 19 +- .../nonstiff/GillFieldStepInterpolatorTest.java | 20 +- .../HighamHall54FieldStepInterpolatorTest.java | 20 +- .../LutherFieldStepInterpolatorTest.java | 18 +- .../MidpointFieldStepInterpolatorTest.java | 20 +- .../ThreeEighthesFieldStepInterpolatorTest.java | 20 +- .../sampling/DummyFieldStepInterpolator.java | 55 +++++ 37 files changed, 981 insertions(+), 568 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java index 07fa9a0..b4f321c 100644 --- a/src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java @@ -307,6 +307,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp } } + AbstractFieldStepInterpolator<T> restricted = interpolator; while (!occurringEvents.isEmpty()) { // handle the chronologically first event @@ -315,11 +316,10 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp iterator.remove(); // get state at event time - final FieldODEStateAndDerivative<T> eventState = interpolator.getInterpolatedState(currentEvent.getEventTime()); + final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime()); // restrict the interpolator to the first part of the step, up to the event - interpolator.setSoftPreviousState(previousState); - interpolator.setSoftCurrentState(eventState); + restricted = restricted.restrictStep(previousState, eventState); // advance all event states to current time for (final FieldEventState<T> state : eventsStates) { @@ -329,7 +329,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp // handle the first part of the step, up to the event for (final FieldStepHandler<T> handler : stepHandlers) { - handler.handleStep(interpolator, isLastStep); + handler.handleStep(restricted, isLastStep); } if (isLastStep) { @@ -351,11 +351,10 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp // prepare handling of the remaining part of the step previousState = eventState; - interpolator.setSoftPreviousState(eventState); - interpolator.setSoftCurrentState(currentState); + restricted = restricted.restrictStep(eventState, currentState); // check if the same event occurs again in the remaining part of the step - if (currentEvent.evaluateStep(interpolator)) { + if (currentEvent.evaluateStep(restricted)) { // the event occurs during the current step occurringEvents.add(currentEvent); } @@ -371,7 +370,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp // handle the remaining part of the step, after all events if any for (FieldStepHandler<T> handler : stepHandlers) { - handler.handleStep(interpolator, isLastStep); + handler.handleStep(restricted, isLastStep); } return currentState; http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java b/src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java index 8f94e51..f26fded 100644 --- a/src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java +++ b/src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java @@ -80,7 +80,7 @@ import org.apache.commons.math4.util.FastMath; */ public class ContinuousOutputFieldModel<T extends RealFieldElement<T>> - implements FieldStepHandler<T> { + implements FieldStepHandler<T> { /** Initial integration time. */ private T initialTime; @@ -156,7 +156,7 @@ public class ContinuousOutputFieldModel<T extends RealFieldElement<T>> } for (FieldStepInterpolator<T> interpolator : model.steps) { - steps.add(interpolator.copy()); + steps.add(interpolator); } index = steps.size() - 1; @@ -201,7 +201,7 @@ public class ContinuousOutputFieldModel<T extends RealFieldElement<T>> forward = interpolator.isForward(); } - steps.add(interpolator.copy()); + steps.add(interpolator); if (isLast) { finalTime = interpolator.getCurrentState().getTime(); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java index a09c1fa..315924a 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; /** @@ -100,8 +101,14 @@ public class ClassicalRungeKuttaFieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected ClassicalRungeKuttaFieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new ClassicalRungeKuttaFieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldEquationsMapper<T> mapper) { + return new ClassicalRungeKuttaFieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java index d096db2..9b70331 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java @@ -62,26 +62,36 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ ClassicalRungeKuttaFieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); - } - - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - ClassicalRungeKuttaFieldStepInterpolator(final ClassicalRungeKuttaFieldStepInterpolator<T> interpolator) { - super(interpolator); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected ClassicalRungeKuttaFieldStepInterpolator<T> doCopy() { - return new ClassicalRungeKuttaFieldStepInterpolator<T>(this); + protected ClassicalRungeKuttaFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new ClassicalRungeKuttaFieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } /** {@inheritDoc} */ @@ -89,9 +99,9 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { - final T one = getField().getOne(); + final T one = time.getField().getOne(); final T oneMinusTheta = one.subtract(theta); final T oneMinus2Theta = one.subtract(theta.multiply(2)); final T coeffDot1 = oneMinusTheta.multiply(oneMinus2Theta); @@ -102,7 +112,7 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>> 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 s = thetaH.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))); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java index e454ca7..b8ec110 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; import org.apache.commons.math4.util.MathUtils; @@ -189,8 +190,13 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected DormandPrince54FieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new DormandPrince54FieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) { + return new DormandPrince54FieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java index 93b51d7..c9495e2 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java @@ -75,11 +75,23 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ DormandPrince54FieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); final T one = field.getOne(); a70 = one.multiply( 35.0).divide( 384.0); a72 = one.multiply( 500.0).divide(1113.0); @@ -94,43 +106,27 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> d6 = one.multiply( 69997945.0).divide( 29380423.0); } - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T> interpolator) { - - super(interpolator); - 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; - - } - /** {@inheritDoc} */ - @Override - protected DormandPrince54FieldStepInterpolator<T> doCopy() { - return new DormandPrince54FieldStepInterpolator<T>(this); + protected DormandPrince54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new DormandPrince54FieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { // interpolate - final T one = getField().getOne(); + final T one = time.getField().getOne(); final T eta = one.subtract(theta); final T twoTheta = theta.multiply(2); final T dot2 = one.subtract(twoTheta); @@ -139,7 +135,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> final T[] interpolatedState; final T[] interpolatedDerivatives; if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - final T f1 = h.multiply(theta); + final T f1 = thetaH; final T f2 = f1.multiply(eta); final T f3 = f2.multiply(theta); final T f4 = f3.multiply(eta); @@ -147,7 +143,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> subtract(f2.multiply(a70.subtract(1))). add(f3.multiply(a70.multiply(2).subtract(1))). add(f4.multiply(d0)); - final T coeff1 = getField().getZero(); + final T coeff1 = time.getField().getZero(); final T coeff2 = f1.multiply(a72). subtract(f2.multiply(a72)). add(f3.multiply(a72.multiply(2))). @@ -169,7 +165,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> subtract(dot2.multiply(a70.subtract(1))). add(dot3.multiply(a70.multiply(2).subtract(1))). add(dot4.multiply(d0)); - final T coeffDot1 = getField().getZero(); + final T coeffDot1 = time.getField().getZero(); final T coeffDot2 = a72. subtract(dot2.multiply(a72)). add(dot3.multiply(a72.multiply(2))). @@ -200,7 +196,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> subtract(f2.multiply(a70.subtract(1))). add(f3.multiply(a70.multiply(2).subtract(1))). add(f4.multiply(d0)); - final T coeff1 = getField().getZero(); + final T coeff1 = time.getField().getZero(); final T coeff2 = f1.multiply(a72). subtract(f2.multiply(a72)). add(f3.multiply(a72.multiply(2))). @@ -222,7 +218,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>> subtract(dot2.multiply(a70.subtract(1))). add(dot3.multiply(a70.multiply(2).subtract(1))). add(dot4.multiply(d0)); - final T coeffDot1 = getField().getZero(); + final T coeffDot1 = time.getField().getZero(); final T coeffDot2 = a72. subtract(dot2.multiply(a72)). add(dot3.multiply(a72.multiply(2))). http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java index 9f8e36c..171b465 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; import org.apache.commons.math4.util.MathUtils; @@ -395,8 +396,13 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected DormandPrince853FieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new DormandPrince853FieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) { + return new DormandPrince853FieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java index 68573fd..9f4b2b0 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java @@ -45,32 +45,43 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ DormandPrince853FieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); - + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); // interpolation weights - d = MathArrays.buildArray(getField(), 7, 16); + d = MathArrays.buildArray(field, 7, 16); // this row is the same as the b array - d[0][ 0] = fraction(104257, 1920240); - d[0][ 1] = getField().getZero(); - d[0][ 2] = getField().getZero(); - d[0][ 3] = getField().getZero(); - d[0][ 4] = 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] = getField().getZero(); - d[0][13] = getField().getZero(); - d[0][14] = getField().getZero(); - d[0][15] = getField().getZero(); + d[0][ 0] = fraction(field, 104257, 1920240); + d[0][ 1] = field.getZero(); + d[0][ 2] = field.getZero(); + d[0][ 3] = field.getZero(); + d[0][ 4] = field.getZero(); + d[0][ 5] = fraction(field, 3399327.0, 763840.0); + d[0][ 6] = fraction(field, 66578432.0, 35198415.0); + d[0][ 7] = fraction(field, -1674902723.0, 288716400.0); + d[0][ 8] = fraction(field, 54980371265625.0, 176692375811392.0); + d[0][ 9] = fraction(field, -734375.0, 4826304.0); + d[0][10] = fraction(field, 171414593.0, 851261400.0); + d[0][11] = fraction(field, 137909.0, 3084480.0); + d[0][12] = field.getZero(); + d[0][13] = field.getZero(); + d[0][14] = field.getZero(); + d[0][15] = field.getZero(); d[1][ 0] = d[0][ 0].negate().add(1); d[1][ 1] = d[0][ 1].negate(); @@ -106,105 +117,97 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>> 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] = getField().getZero(); - d[3][ 2] = getField().getZero(); - d[3][ 3] = getField().getZero(); - d[3][ 4] = 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[3][ 0] = fraction(field, -17751989329.0, 2106076560.0); + d[3][ 1] = field.getZero(); + d[3][ 2] = field.getZero(); + d[3][ 3] = field.getZero(); + d[3][ 4] = field.getZero(); + d[3][ 5] = fraction(field, 4272954039.0, 7539864640.0); + d[3][ 6] = fraction(field, -118476319744.0, 38604839385.0); + d[3][ 7] = fraction(field, 755123450731.0, 316657731600.0); + d[3][ 8] = fraction(field, 3692384461234828125.0, 1744130441634250432.0); + d[3][ 9] = fraction(field, -4612609375.0, 5293382976.0); + d[3][10] = fraction(field, 2091772278379.0, 933644586600.0); + d[3][11] = fraction(field, 2136624137.0, 3382989120.0); + d[3][12] = fraction(field, -126493.0, 1421424.0); + d[3][13] = fraction(field, 98350000.0, 5419179.0); + d[3][14] = fraction(field, -18878125.0, 2053168.0); + d[3][15] = fraction(field, -1944542619.0, 438351368.0); - d[4][ 0] = fraction( 32941697297.0, 3159114840.0); - d[4][ 1] = getField().getZero(); - d[4][ 2] = getField().getZero(); - d[4][ 3] = getField().getZero(); - d[4][ 4] = 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[4][ 0] = fraction(field, 32941697297.0, 3159114840.0); + d[4][ 1] = field.getZero(); + d[4][ 2] = field.getZero(); + d[4][ 3] = field.getZero(); + d[4][ 4] = field.getZero(); + d[4][ 5] = fraction(field, 456696183123.0, 1884966160.0); + d[4][ 6] = fraction(field, 19132610714624.0, 115814518155.0); + d[4][ 7] = fraction(field, -177904688592943.0, 474986597400.0); + d[4][ 8] = fraction(field, -4821139941836765625.0, 218016305204281304.0); + d[4][ 9] = fraction(field, 30702015625.0, 3970037232.0); + d[4][10] = fraction(field, -85916079474274.0, 2800933759800.0); + d[4][11] = fraction(field, -5919468007.0, 634310460.0); + d[4][12] = fraction(field, 2479159.0, 157936.0); + d[4][13] = fraction(field, -18750000.0, 602131.0); + d[4][14] = fraction(field, -19203125.0, 2053168.0); + d[4][15] = fraction(field, 15700361463.0, 438351368.0); - d[5][ 0] = fraction( 12627015655.0, 631822968.0); - d[5][ 1] = getField().getZero(); - d[5][ 2] = getField().getZero(); - d[5][ 3] = getField().getZero(); - d[5][ 4] = 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[5][ 0] = fraction(field, 12627015655.0, 631822968.0); + d[5][ 1] = field.getZero(); + d[5][ 2] = field.getZero(); + d[5][ 3] = field.getZero(); + d[5][ 4] = field.getZero(); + d[5][ 5] = fraction(field, -72955222965.0, 188496616.0); + d[5][ 6] = fraction(field, -13145744952320.0, 69488710893.0); + d[5][ 7] = fraction(field, 30084216194513.0, 56998391688.0); + d[5][ 8] = fraction(field, -296858761006640625.0, 25648977082856624.0); + d[5][ 9] = fraction(field, 569140625.0, 82709109.0); + d[5][10] = fraction(field, -18684190637.0, 18672891732.0); + d[5][11] = fraction(field, 69644045.0, 89549712.0); + d[5][12] = fraction(field, -11847025.0, 4264272.0); + d[5][13] = fraction(field, -978650000.0, 16257537.0); + d[5][14] = fraction(field, 519371875.0, 6159504.0); + d[5][15] = fraction(field, 5256837225.0, 438351368.0); - d[6][ 0] = fraction( -450944925.0, 17550638.0); - d[6][ 1] = getField().getZero(); - d[6][ 2] = getField().getZero(); - d[6][ 3] = getField().getZero(); - d[6][ 4] = 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); + d[6][ 0] = fraction(field, -450944925.0, 17550638.0); + d[6][ 1] = field.getZero(); + d[6][ 2] = field.getZero(); + d[6][ 3] = field.getZero(); + d[6][ 4] = field.getZero(); + d[6][ 5] = fraction(field, -14532122925.0, 94248308.0); + d[6][ 6] = fraction(field, -595876966400.0, 2573655959.0); + d[6][ 7] = fraction(field, 188748653015.0, 527762886.0); + d[6][ 8] = fraction(field, 2545485458115234375.0, 27252038150535163.0); + d[6][ 9] = fraction(field, -1376953125.0, 36759604.0); + d[6][10] = fraction(field, 53995596795.0, 518691437.0); + d[6][11] = fraction(field, 210311225.0, 7047894.0); + d[6][12] = fraction(field, -1718875.0, 39484.0); + d[6][13] = fraction(field, 58000000.0, 602131.0); + d[6][14] = fraction(field, -1546875.0, 39484.0); + d[6][15] = fraction(field, -1262172375.0, 8429834.0); } - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - DormandPrince853FieldStepInterpolator(final DormandPrince853FieldStepInterpolator<T> interpolator) { - - super(interpolator); - - d = MathArrays.buildArray(getField(), 4, -1); - for (int i = 0; i < d.length; ++i) { - d[i] = interpolator.d[i].clone(); - } - + /** {@inheritDoc} */ + protected DormandPrince853FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new DormandPrince853FieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } /** Create a fraction. + * @param field field to which the elements belong * @param p numerator * @param q denominator * @return p/q computed in the instance field */ - private T fraction(final double p, final double q) { - return getField().getOne().multiply(p).divide(q); - } - - /** {@inheritDoc} */ - @Override - protected DormandPrince853FieldStepInterpolator<T> doCopy() { - return new DormandPrince853FieldStepInterpolator<T>(this); + private T fraction(final Field<T> field, final double p, final double q) { + return field.getZero().add(p).divide(q); } /** {@inheritDoc} */ @@ -212,10 +215,10 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>> @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) + final T thetaH, final T oneMinusThetaH) throws MaxCountExceededException { - final T one = getField().getOne(); + final T one = time.getField().getOne(); final T eta = one.subtract(theta); final T twoTheta = theta.multiply(2); final T theta2 = theta.multiply(theta); @@ -230,15 +233,15 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>> if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - final T f0 = theta.multiply(h); + final T f0 = thetaH; 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); - final T[] p = MathArrays.buildArray(getField(), 16); - final T[] q = MathArrays.buildArray(getField(), 16); + final T[] p = MathArrays.buildArray(time.getField(), 16); + final T[] q = MathArrays.buildArray(time.getField(), 16); for (int i = 0; i < p.length; ++i) { p[i] = f0.multiply(d[0][i]). add(f1.multiply(d[1][i])). @@ -267,8 +270,8 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>> final T f4 = f3.multiply(theta); final T f5 = f4.multiply(eta); final T f6 = f5.multiply(theta); - final T[] p = MathArrays.buildArray(getField(), 16); - final T[] q = MathArrays.buildArray(getField(), 16); + final T[] p = MathArrays.buildArray(time.getField(), 16); + final T[] q = MathArrays.buildArray(time.getField(), 16); for (int i = 0; i < p.length; ++i) { p[i] = f0.multiply(d[0][i]). add(f1.multiply(d[1][i])). http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java index 13bb070..c6e2444 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java @@ -183,10 +183,15 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme /** Create an interpolator. * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step * @param mapper equations mapper for the all equations * @return external weights for the high order method from Butcher array */ - protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, + protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, FieldEquationsMapper<T> mapper); /** Get the order of the method. * @return order of the method @@ -226,11 +231,6 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1); final T[] yTmp = MathArrays.buildArray(getField(), y0.length); - // set up an interpolator sharing the integrator arrays - final RungeKuttaFieldStepInterpolator<T> interpolator = - createInterpolator(forward, equations.getMapper()); - interpolator.storeState(stepStart); - // set up integration control objects T hNew = getField().getZero(); boolean firstTime = true; @@ -239,8 +239,6 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme isLastStep = false; do { - interpolator.shift(); - // iterate over step size, ensuring local normalized error is smaller than 1 T error = getField().getZero().add(10); while (error.subtract(1.0).getReal() >= 0) { @@ -314,16 +312,12 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme 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); + stepStart = acceptStep(createInterpolator(forward, yDotK, stepStart, stateTmp, equations.getMapper()), + finalTime); if (!isLastStep) { - // prepare next step - interpolator.storeState(stepStart); - // stepsize control for next step final T factor = MathUtils.min(maxGrowth, MathUtils.max(minReduction, safety.multiply(error.pow(exp)))); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java index 96403e5..ce87431 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; /** @@ -85,8 +86,14 @@ public class EulerFieldIntegrator<T extends RealFieldElement<T>> extends RungeKu /** {@inheritDoc} */ @Override protected EulerFieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new EulerFieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldEquationsMapper<T> mapper) { + return new EulerFieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java index 9046118..133e540 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java @@ -52,26 +52,36 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ EulerFieldStepInterpolator(final Field<T> field, final boolean forward, - final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); - } - - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - EulerFieldStepInterpolator(final EulerFieldStepInterpolator<T> interpolator) { - super(interpolator); + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, + final FieldEquationsMapper<T> mapper) { + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected EulerFieldStepInterpolator<T> doCopy() { - return new EulerFieldStepInterpolator<T>(this); + protected EulerFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new EulerFieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } /** {@inheritDoc} */ @@ -79,15 +89,15 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>> @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { final T[] interpolatedState; final T[] interpolatedDerivatives; if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - interpolatedState = previousStateLinearCombination(theta.multiply(h)); - interpolatedDerivatives = derivativeLinearCombination(getField().getOne()); + interpolatedState = previousStateLinearCombination(thetaH); + interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne()); } else { interpolatedState = currentStateLinearCombination(oneMinusThetaH.negate()); - interpolatedDerivatives = derivativeLinearCombination(getField().getOne()); + interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne()); } return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java index 6996d05..1fc6688 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; @@ -110,8 +111,14 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected GillFieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new GillFieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldEquationsMapper<T> mapper) { + return new GillFieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java index 62590e7..fdb16cb 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java @@ -67,42 +67,49 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ GillFieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); final T sqrt = field.getZero().add(0.5).sqrt(); one_minus_inv_sqrt_2 = field.getOne().subtract(sqrt); one_plus_inv_sqrt_2 = field.getOne().add(sqrt); } - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - GillFieldStepInterpolator(final GillFieldStepInterpolator<T> interpolator) { - super(interpolator); - one_minus_inv_sqrt_2 = interpolator.one_minus_inv_sqrt_2; - one_plus_inv_sqrt_2 = interpolator.one_plus_inv_sqrt_2; - } - /** {@inheritDoc} */ - @Override - protected GillFieldStepInterpolator<T> doCopy() { - return new GillFieldStepInterpolator<T>(this); + protected GillFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new GillFieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { - final T one = getField().getOne(); + final T one = time.getField().getOne(); final T twoTheta = theta.multiply(2); final T fourTheta2 = twoTheta.multiply(twoTheta); final T coeffDot1 = theta.multiply(twoTheta.subtract(3)).add(1); @@ -114,7 +121,7 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>> final T[] interpolatedDerivatives; if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - final T s = theta.multiply(h).divide(6.0); + final T s = thetaH.divide(6.0); final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2)); final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6)); final T coeff2 = c23.multiply(one_minus_inv_sqrt_2); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java index 0e0c3b8..0f8353a 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; import org.apache.commons.math4.util.MathUtils; @@ -164,8 +165,13 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected HighamHall54FieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new HighamHall54FieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) { + return new HighamHall54FieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java index c278043..c5402c2 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java @@ -38,38 +38,47 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations - */ + */ HighamHall54FieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); - } - - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - HighamHall54FieldStepInterpolator(final HighamHall54FieldStepInterpolator<T> interpolator) { - super(interpolator); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected HighamHall54FieldStepInterpolator<T> doCopy() { - return new HighamHall54FieldStepInterpolator<T>(this); + protected HighamHall54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new HighamHall54FieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, 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 = getField().getZero(); + 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)); @@ -78,19 +87,19 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>> final T[] interpolatedDerivatives; 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 = 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))); + final T b0 = thetaH.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 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0))); + final T b3 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 ))); + final T b4 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0))); + final T b5 = thetaH.multiply(theta.multiply(theta.multiply( 5.0 / 12.0 ).add( -5.0 / 16.0))); 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 h = thetaH.divide(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 = getField().getZero(); + 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)); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java index b8d78d9..1c3301c 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; @@ -135,8 +136,14 @@ public class LutherFieldIntegrator<T extends RealFieldElement<T>> /** {@inheritDoc} */ @Override protected LutherFieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new LutherFieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldEquationsMapper<T> mapper) { + return new LutherFieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java index a5dc565..22bc403 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java @@ -83,11 +83,23 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ LutherFieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); final T q = field.getZero().add(21).sqrt(); c5a = q.multiply( -49).add( -49); c5b = q.multiply( 287).add( 392); @@ -103,45 +115,27 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> d6a = q.multiply( -49).add( 49); d6b = q.multiply( 847).add(-1372); d6c = q.multiply(-1029).add( 2254); - - } - - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - LutherFieldStepInterpolator(final LutherFieldStepInterpolator<T> interpolator) { - super(interpolator); - c5a = interpolator.c5a; - c5b = interpolator.c5b; - c5c = interpolator.c5c; - c5d = interpolator.c5d; - c6a = interpolator.c6a; - c6b = interpolator.c6b; - c6c = interpolator.c6c; - c6d = interpolator.c6d; - d5a = interpolator.d5a; - d5b = interpolator.d5b; - d5c = interpolator.d5c; - d6a = interpolator.d6a; - d6b = interpolator.d6b; - d6c = interpolator.d6c; } /** {@inheritDoc} */ - @Override - protected LutherFieldStepInterpolator<T> doCopy() { - return new LutherFieldStepInterpolator<T>(this); + protected LutherFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new LutherFieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { // the coefficients below have been computed by solving the // order conditions from a theorem from Butcher (1963), using @@ -187,7 +181,7 @@ 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); - final T coeffDot2 = getField().getZero(); + final T coeffDot2 = time.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))); @@ -198,9 +192,9 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - final T s = theta.multiply(h); + final T s = thetaH; final T coeff1 = s.multiply(theta.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 = getField().getZero(); + final T coeff2 = time.getField().getZero(); final T coeff3 = s.multiply(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 = s.multiply(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 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300)))); @@ -212,7 +206,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>> final T s = oneMinusThetaH; final T coeff1 = s.multiply(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)); - final T coeff2 = getField().getZero(); + final T coeff2 = time.getField().getZero(); final T coeff3 = s.multiply(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 = s.multiply(theta.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(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0)); http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java index 2a2e5e6..af5a867 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java @@ -20,6 +20,7 @@ package org.apache.commons.math4.ode.nonstiff; import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.util.MathArrays; /** @@ -85,8 +86,14 @@ public class MidpointFieldIntegrator<T extends RealFieldElement<T>> extends Rung /** {@inheritDoc} */ @Override protected MidpointFieldStepInterpolator<T> - createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) { - return new MidpointFieldStepInterpolator<T>(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldEquationsMapper<T> mapper) { + return new MidpointFieldStepInterpolator<T>(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } http://git-wip-us.apache.org/repos/asf/commons-math/blob/771eb6a6/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java index 23e00b6..5a01a35 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java @@ -54,44 +54,53 @@ class MidpointFieldStepInterpolator<T extends RealFieldElement<T>> /** Simple constructor. * @param field field to which the time and state vector elements belong * @param forward integration direction indicator + * @param yDotK slopes at the intermediate points + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param mapper equations mapper for the all equations */ - MidpointFieldStepInterpolator(Field<T> field, final boolean forward, - final FieldEquationsMapper<T> mapper) { - super(field, forward, mapper); - } - - /** Copy constructor. - * @param interpolator interpolator to copy from. The copy is a deep - * copy: its arrays are separated from the original arrays of the - * instance - */ - MidpointFieldStepInterpolator(final MidpointFieldStepInterpolator<T> interpolator) { - super(interpolator); + MidpointFieldStepInterpolator(final Field<T> field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, + final FieldEquationsMapper<T> mapper) { + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected MidpointFieldStepInterpolator<T> doCopy() { - return new MidpointFieldStepInterpolator<T>(this); + protected MidpointFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative<T> newGlobalPreviousState, + final FieldODEStateAndDerivative<T> newGlobalCurrentState, + final FieldODEStateAndDerivative<T> newSoftPreviousState, + final FieldODEStateAndDerivative<T> newSoftCurrentState, + final FieldEquationsMapper<T> newMapper) { + return new MidpointFieldStepInterpolator<T>(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { final T coeffDot2 = theta.multiply(2); - final T coeffDot1 = getField().getOne().subtract(coeffDot2); + final T coeffDot1 = time.getField().getOne().subtract(coeffDot2); final T[] interpolatedState; final T[] interpolatedDerivatives; if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { final T coeff1 = theta.multiply(oneMinusThetaH); - final T coeff2 = theta.multiply(theta).multiply(h); + final T coeff2 = theta.multiply(thetaH); interpolatedState = previousStateLinearCombination(coeff1, coeff2); interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2); } else {
