From 771eb6a606491375a61fee1aa537025edafacc34 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 6 Jan 2016 13:22:39 +0100 Subject: [PATCH] Use immutable step interpolators. --- .../math4/ode/AbstractFieldIntegrator.java | 15 +- .../math4/ode/ContinuousOutputFieldModel.java | 6 +- .../ClassicalRungeKuttaFieldIntegrator.java | 11 +- ...ssicalRungeKuttaFieldStepInterpolator.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 ++-- .../ode/nonstiff/GillFieldIntegrator.java | 11 +- .../nonstiff/GillFieldStepInterpolator.java | 45 ++-- .../nonstiff/HighamHall54FieldIntegrator.java | 10 +- .../HighamHall54FieldStepInterpolator.java | 61 +++-- .../ode/nonstiff/LutherFieldIntegrator.java | 11 +- .../nonstiff/LutherFieldStepInterpolator.java | 62 +++-- .../ode/nonstiff/MidpointFieldIntegrator.java | 11 +- .../MidpointFieldStepInterpolator.java | 47 ++-- .../nonstiff/RungeKuttaFieldIntegrator.java | 22 +- .../RungeKuttaFieldStepInterpolator.java | 83 ++++--- .../ThreeEighthesFieldIntegrator.java | 11 +- .../ThreeEighthesFieldStepInterpolator.java | 41 ++-- .../AbstractFieldStepInterpolator.java | 145 ++++------- .../ode/sampling/FieldStepInterpolator.java | 11 - .../ode/ContinuousOutputFieldModelTest.java | 219 +++++++++++++++++ ...ctRungeKuttaFieldStepInterpolatorTest.java | 71 ++---- ...calRungKuttaFieldStepInterpolatorTest.java | 20 +- ...mandPrince54FieldStepInterpolatorTest.java | 20 +- ...andPrince853FieldStepInterpolatorTest.java | 18 +- .../EulerFieldStepInterpolatorTest.java | 19 +- .../GillFieldStepInterpolatorTest.java | 20 +- ...HighamHall54FieldStepInterpolatorTest.java | 20 +- .../LutherFieldStepInterpolatorTest.java | 18 +- .../MidpointFieldStepInterpolatorTest.java | 20 +- ...hreeEighthesFieldStepInterpolatorTest.java | 20 +- .../sampling/DummyFieldStepInterpolator.java | 55 +++++ 37 files changed, 984 insertions(+), 571 deletions(-) create mode 100644 src/test/java/org/apache/commons/math4/ode/ContinuousOutputFieldModelTest.java create mode 100644 src/test/java/org/apache/commons/math4/ode/sampling/DummyFieldStepInterpolator.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 07fa9a050..b4f321ce7 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> imp } } + AbstractFieldStepInterpolator restricted = interpolator; while (!occurringEvents.isEmpty()) { // handle the chronologically first event @@ -315,11 +316,10 @@ public abstract class AbstractFieldIntegrator> imp iterator.remove(); // get state at event time - final FieldODEStateAndDerivative eventState = interpolator.getInterpolatedState(currentEvent.getEventTime()); + final FieldODEStateAndDerivative 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 state : eventsStates) { @@ -329,7 +329,7 @@ public abstract class AbstractFieldIntegrator> imp // handle the first part of the step, up to the event for (final FieldStepHandler handler : stepHandlers) { - handler.handleStep(interpolator, isLastStep); + handler.handleStep(restricted, isLastStep); } if (isLastStep) { @@ -351,11 +351,10 @@ public abstract class AbstractFieldIntegrator> 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> imp // handle the remaining part of the step, after all events if any for (FieldStepHandler handler : stepHandlers) { - handler.handleStep(interpolator, isLastStep); + handler.handleStep(restricted, isLastStep); } return currentState; 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 8f94e5178..f26fded7a 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> - implements FieldStepHandler { + implements FieldStepHandler { /** Initial integration time. */ private T initialTime; @@ -156,7 +156,7 @@ public class ContinuousOutputFieldModel> } for (FieldStepInterpolator interpolator : model.steps) { - steps.add(interpolator.copy()); + steps.add(interpolator); } index = steps.size() - 1; @@ -201,7 +201,7 @@ public class ContinuousOutputFieldModel> forward = interpolator.isForward(); } - steps.add(interpolator.copy()); + steps.add(interpolator); if (isLast) { finalTime = interpolator.getCurrentState().getTime(); 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 a09c1fa9e..315924a90 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> /** {@inheritDoc} */ @Override protected ClassicalRungeKuttaFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new ClassicalRungeKuttaFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new ClassicalRungeKuttaFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } 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 d096db21f..9b7033130 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> /** 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 field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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 interpolator) { - super(interpolator); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected ClassicalRungeKuttaFieldStepInterpolator doCopy() { - return new ClassicalRungeKuttaFieldStepInterpolator(this); + protected ClassicalRungeKuttaFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new ClassicalRungeKuttaFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } /** {@inheritDoc} */ @@ -89,9 +99,9 @@ class ClassicalRungeKuttaFieldStepInterpolator> @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> 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))); 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 e454ca75a..b8ec110e8 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> /** {@inheritDoc} */ @Override protected DormandPrince54FieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new DormandPrince54FieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, final FieldEquationsMapper mapper) { + return new DormandPrince54FieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ 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 93b51d769..c9495e2c2 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> /** 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 field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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> 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 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 doCopy() { - return new DormandPrince54FieldStepInterpolator(this); + protected DormandPrince54FieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new DormandPrince54FieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> 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> 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> 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> 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> 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))). 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 9f8e36c7a..171b46524 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> /** {@inheritDoc} */ @Override protected DormandPrince853FieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new DormandPrince853FieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, final FieldEquationsMapper mapper) { + return new DormandPrince853FieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ 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 68573fd76..9f4b2b07d 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> /** 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 field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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> 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 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 create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new DormandPrince853FieldStepInterpolator(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 doCopy() { - return new DormandPrince853FieldStepInterpolator(this); + private T fraction(final Field field, final double p, final double q) { + return field.getZero().add(p).divide(q); } /** {@inheritDoc} */ @@ -212,10 +215,10 @@ class DormandPrince853FieldStepInterpolator> @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> 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> 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])). 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 13bb0704e..c6e244461 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 createInterpolator(boolean forward, + protected abstract RungeKuttaFieldStepInterpolator createInterpolator(boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, FieldEquationsMapper mapper); /** Get the order of the method. * @return order of the method @@ -226,11 +231,6 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator 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= 0) { @@ -314,16 +312,12 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator stateTmp = new FieldODEStateAndDerivative(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)))); 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 96403e589..ce8743111 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> extends RungeKu /** {@inheritDoc} */ @Override protected EulerFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new EulerFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new EulerFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } 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 90461189c..133e5405c 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> /** 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 field, final boolean forward, - final FieldEquationsMapper 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 interpolator) { - super(interpolator); + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, + final FieldEquationsMapper mapper) { + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected EulerFieldStepInterpolator doCopy() { - return new EulerFieldStepInterpolator(this); + protected EulerFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new EulerFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } /** {@inheritDoc} */ @@ -79,15 +89,15 @@ class EulerFieldStepInterpolator> @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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(time, interpolatedState, interpolatedDerivatives); 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 6996d0582..1fc66887c 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> /** {@inheritDoc} */ @Override protected GillFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new GillFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new GillFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } 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 62590e75a..fdb16cbff 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> /** 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 field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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 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 doCopy() { - return new GillFieldStepInterpolator(this); + protected GillFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new GillFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> 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); 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 0e0c3b8b9..0f8353af4 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> /** {@inheritDoc} */ @Override protected HighamHall54FieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new HighamHall54FieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, final FieldEquationsMapper mapper) { + return new HighamHall54FieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } /** {@inheritDoc} */ 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 c27804314..c5402c202 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> /** 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 field, final boolean forward, - final FieldEquationsMapper 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 interpolator) { - super(interpolator); + HighamHall54FieldStepInterpolator(final Field field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, + final FieldEquationsMapper mapper) { + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected HighamHall54FieldStepInterpolator doCopy() { - return new HighamHall54FieldStepInterpolator(this); + protected HighamHall54FieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new HighamHall54FieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> 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)); 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 b8d78d91e..1c3301c8b 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> /** {@inheritDoc} */ @Override protected LutherFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new LutherFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new LutherFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } 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 a5dc565d9..22bc403f2 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> /** 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 field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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> 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 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 doCopy() { - return new LutherFieldStepInterpolator(this); + protected LutherFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new LutherFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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> // 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> 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> 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)); 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 2a2e5e6d6..af5a867f3 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> extends Rung /** {@inheritDoc} */ @Override protected MidpointFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new MidpointFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new MidpointFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } 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 23e00b6c3..5a01a352f 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> /** 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 field, final boolean forward, - final FieldEquationsMapper 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 interpolator) { - super(interpolator); + MidpointFieldStepInterpolator(final Field field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, + final FieldEquationsMapper mapper) { + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected MidpointFieldStepInterpolator doCopy() { - return new MidpointFieldStepInterpolator(this); + protected MidpointFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new MidpointFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper 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 { diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java index 3a5f9aacb..d052f4b29 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java @@ -99,10 +99,15 @@ public abstract class RungeKuttaFieldIntegrator> /** 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 createInterpolator(boolean forward, + protected abstract RungeKuttaFieldStepInterpolator createInterpolator(boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, FieldEquationsMapper mapper); /** {@inheritDoc} */ @@ -124,11 +129,6 @@ public abstract class RungeKuttaFieldIntegrator> 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 interpolator = - createInterpolator(forward, equations.getMapper()); - interpolator.storeState(stepStart); - // set up integration control objects if (forward) { if (stepStart.getTime().add(step).subtract(finalTime).getReal() >= 0) { @@ -148,8 +148,6 @@ public abstract class RungeKuttaFieldIntegrator> isLastStep = false; do { - interpolator.shift(); - // first stage y = equations.getMapper().mapState(stepStart); yDotK[0] = equations.getMapper().mapDerivative(stepStart); @@ -182,16 +180,12 @@ public abstract class RungeKuttaFieldIntegrator> final FieldODEStateAndDerivative stateTmp = new FieldODEStateAndDerivative(stepEnd, yTmp, yDotTmp); // discrete events handling - 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 nextT = stepStart.getTime().add(stepSize); final boolean nextIsLast = forward ? diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.java index 26d8fb990..9595ee247 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.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.ode.sampling.AbstractFieldStepInterpolator; import org.apache.commons.math4.util.MathArrays; @@ -40,57 +41,63 @@ abstract class RungeKuttaFieldStepInterpolator> private final Field field; /** Slopes at the intermediate points. */ - private T[][] yDotK; + private final T[][] yDotK; /** 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 */ protected RungeKuttaFieldStepInterpolator(final Field field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper mapper) { - super(forward, mapper); + super(forward, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper); this.field = field; - this.yDotK = null; - } - - /** Copy constructor. - *

The copy is a deep copy: its arrays are separated from the - * original arrays of the instance.

- - * @param interpolator interpolator to copy from. - - */ - RungeKuttaFieldStepInterpolator(final RungeKuttaFieldStepInterpolator interpolator) { - - super(interpolator); - field = interpolator.field; - - if (yDotK != null) { - yDotK = MathArrays.buildArray(field, interpolator.yDotK.length, -1); - for (int k = 0; k < yDotK.length; ++k) { - yDotK[k] = interpolator.yDotK[k].clone(); - } - - } else { - yDotK = null; + this.yDotK = MathArrays.buildArray(field, yDotK.length, -1); + for (int i = 0; i < yDotK.length; ++i) { + this.yDotK[i] = yDotK[i].clone(); } - } - /** Get the field to which the time and state vector elements belong. - * @return to which the time and state vector elements belong + /** {@inheritDoc} */ + protected RungeKuttaFieldStepInterpolator create(boolean newForward, + FieldODEStateAndDerivative newGlobalPreviousState, + FieldODEStateAndDerivative newGlobalCurrentState, + FieldODEStateAndDerivative newSoftPreviousState, + FieldODEStateAndDerivative newSoftCurrentState, + FieldEquationsMapper newMapper) { + return create(field, newForward, yDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); + } + + /** Create a new instance. + * @param newField field to which the time and state vector elements belong + * @param newForward integration direction indicator + * @param newYDotK slopes at the intermediate points + * @param newGlobalPreviousState start of the global step + * @param newGlobalCurrentState end of the global step + * @param newSoftPreviousState start of the restricted step + * @param newSoftCurrentState end of the restricted step + * @param newMapper equations mapper for the all equations + * @return a new instance */ - protected Field getField() { - return field; - } - - /** Store the slopes at the intermediate points. - * @param slopes slopes at the intermediate points - */ - void setSlopes(final T[][] slopes) { - this.yDotK = slopes.clone(); - } + protected abstract RungeKuttaFieldStepInterpolator create(Field newField, boolean newForward, T[][] newYDotK, + FieldODEStateAndDerivative newGlobalPreviousState, + FieldODEStateAndDerivative newGlobalCurrentState, + FieldODEStateAndDerivative newSoftPreviousState, + FieldODEStateAndDerivative newSoftCurrentState, + FieldEquationsMapper newMapper); /** Compute a state by linear combination added to previous state. * @param coefficients coefficients to apply to the method staged derivatives diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java index b82629061..255b2dafa 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.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; /** @@ -99,8 +100,14 @@ public class ThreeEighthesFieldIntegrator> /** {@inheritDoc} */ @Override protected ThreeEighthesFieldStepInterpolator - createInterpolator(final boolean forward, final FieldEquationsMapper mapper) { - return new ThreeEighthesFieldStepInterpolator(getField(), forward, mapper); + createInterpolator(final boolean forward, T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldEquationsMapper mapper) { + return new ThreeEighthesFieldStepInterpolator(getField(), forward, yDotK, + globalPreviousState, globalCurrentState, + globalPreviousState, globalCurrentState, + mapper); } } diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java index bf0ad4621..6faffd1cc 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java @@ -64,35 +64,44 @@ class ThreeEighthesFieldStepInterpolator> /** 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 */ ThreeEighthesFieldStepInterpolator(final Field field, final boolean forward, + final T[][] yDotK, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper 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 - */ - ThreeEighthesFieldStepInterpolator(final ThreeEighthesFieldStepInterpolator interpolator) { - super(interpolator); + super(field, forward, yDotK, + globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, + mapper); } /** {@inheritDoc} */ - @Override - protected ThreeEighthesFieldStepInterpolator doCopy() { - return new ThreeEighthesFieldStepInterpolator(this); + protected ThreeEighthesFieldStepInterpolator create(final Field newField, final boolean newForward, final T[][] newYDotK, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new ThreeEighthesFieldStepInterpolator(newField, newForward, newYDotK, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); } - /** {@inheritDoc} */ @SuppressWarnings("unchecked") @Override protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(final FieldEquationsMapper mapper, final T time, final T theta, - final T oneMinusThetaH) { + final T thetaH, final T oneMinusThetaH) { final T coeffDot3 = theta.multiply(0.75); final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1); @@ -102,7 +111,7 @@ class ThreeEighthesFieldStepInterpolator> final T[] interpolatedDerivatives; if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { - final T s = theta.multiply(h).divide(8); + final T s = thetaH.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); diff --git a/src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java index 00740f438..0486cefd8 100644 --- a/src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java @@ -40,118 +40,76 @@ import org.apache.commons.math4.ode.FieldODEStateAndDerivative; public abstract class AbstractFieldStepInterpolator> implements FieldStepInterpolator { - /** Current time step. */ - protected T h; - /** Global previous state. */ - private FieldODEStateAndDerivative globalPreviousState; + private final FieldODEStateAndDerivative globalPreviousState; /** Global current state. */ - private FieldODEStateAndDerivative globalCurrentState; + private final FieldODEStateAndDerivative globalCurrentState; /** Soft previous state. */ - private FieldODEStateAndDerivative softPreviousState; + private final FieldODEStateAndDerivative softPreviousState; /** Soft current state. */ - private FieldODEStateAndDerivative softCurrentState; + private final FieldODEStateAndDerivative softCurrentState; /** integration direction. */ - private boolean forward; + private final boolean forward; /** Mapper for ODE equations primary and secondary components. */ private FieldEquationsMapper mapper; /** Simple constructor. * @param isForward integration direction indicator + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step * @param equationsMapper mapper for ODE equations primary and secondary components */ protected AbstractFieldStepInterpolator(final boolean isForward, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, final FieldEquationsMapper equationsMapper) { - globalPreviousState = null; - globalCurrentState = null; - softPreviousState = null; - softCurrentState = null; - h = null; - this.forward = isForward; - this.mapper = equationsMapper; + this.forward = isForward; + this.globalPreviousState = globalPreviousState; + this.globalCurrentState = globalCurrentState; + this.softPreviousState = softPreviousState; + this.softCurrentState = softCurrentState; + this.mapper = equationsMapper; } - /** Copy constructor. - *

The copy is a deep copy: its arrays are separated from the - * original arrays of the instance.

- * @param interpolator interpolator to copy from. - */ - protected AbstractFieldStepInterpolator(final AbstractFieldStepInterpolator interpolator) { - - globalPreviousState = interpolator.globalPreviousState; - globalCurrentState = interpolator.globalCurrentState; - softPreviousState = interpolator.softPreviousState; - softCurrentState = interpolator.softCurrentState; - h = interpolator.h; - forward = interpolator.forward; - mapper = interpolator.mapper; - - } - - /** {@inheritDoc} */ - public FieldStepInterpolator copy() throws MaxCountExceededException { - - // create the new independent instance - return doCopy(); - - } - - /** Really copy the instance. - * @return a copy of the instance - */ - protected abstract FieldStepInterpolator doCopy(); - - /** Shift one step forward. - * Copy the current time into the previous time, hence preparing the - * interpolator for future calls to {@link #storeTime storeTime} - */ - public void shift() { - globalPreviousState = globalCurrentState; - softPreviousState = globalPreviousState; - softCurrentState = globalCurrentState; - } - - /** Store the current step state. - * @param state current state - */ - public void storeState(final FieldODEStateAndDerivative state) { - globalCurrentState = state; - softCurrentState = globalCurrentState; - if (globalPreviousState != null) { - h = globalCurrentState.getTime().subtract(globalPreviousState.getTime()); - } - } - - /** Restrict step range to a limited part of the global step. + /** Create a new restricted version of the instance. *

- * This method can be used to restrict a step and make it appear - * as if the original step was smaller. Calling this method - * only changes the value returned by {@link #getPreviousState()}, - * it does not change any other property + * The instance is not changed at all. *

- * @param softPreviousState start of the restricted step + * @param previousState start of the restricted step + * @param currentState end of the restricted step + * @return restricted version of the instance + * @see #getPreviousState() + * @see #getCurrentState() */ - public void setSoftPreviousState(final FieldODEStateAndDerivative softPreviousState) { - this.softPreviousState = softPreviousState; + public AbstractFieldStepInterpolator restrictStep(final FieldODEStateAndDerivative previousState, + final FieldODEStateAndDerivative currentState) { + return create(forward, globalPreviousState, globalCurrentState, previousState, currentState, mapper); } - /** Restrict step range to a limited part of the global step. - *

- * This method can be used to restrict a step and make it appear - * as if the original step was smaller. Calling this method - * only changes the value returned by {@link #getCurrentState()}, - * it does not change any other property - *

- * @param softCurrentState end of the restricted step + /** Create a new instance. + * @param newForward integration direction indicator + * @param newGlobalPreviousState start of the global step + * @param newGlobalCurrentState end of the global step + * @param newSoftPreviousState start of the restricted step + * @param newSoftCurrentState end of the restricted step + * @param newMapper equations mapper for the all equations + * @return a new instance */ - public void setSoftCurrentState(final FieldODEStateAndDerivative softCurrentState) { - this.softCurrentState = softCurrentState; - } + protected abstract AbstractFieldStepInterpolator create(boolean newForward, + FieldODEStateAndDerivative newGlobalPreviousState, + FieldODEStateAndDerivative newGlobalCurrentState, + FieldODEStateAndDerivative newSoftPreviousState, + FieldODEStateAndDerivative newSoftCurrentState, + FieldEquationsMapper newMapper); /** * Get the previous global grid point state. @@ -169,25 +127,22 @@ public abstract class AbstractFieldStepInterpolator getPreviousState() { return softPreviousState; } - /** {@inheritDoc} - * @see #setSoftCurrentState(FieldODEStateAndDerivative) - */ + /** {@inheritDoc} */ public FieldODEStateAndDerivative getCurrentState() { return softCurrentState; } /** {@inheritDoc} */ public FieldODEStateAndDerivative getInterpolatedState(final T time) { + final T thetaH = time.subtract(globalPreviousState.getTime()); final T oneMinusThetaH = globalCurrentState.getTime().subtract(time); - final T theta = (h.getReal() == 0) ? h.getField().getZero() : h.subtract(oneMinusThetaH).divide(h); - return computeInterpolatedStateAndDerivatives(mapper, time, theta, oneMinusThetaH); + final T theta = thetaH.divide(globalCurrentState.getTime().subtract(globalPreviousState.getTime())); + return computeInterpolatedStateAndDerivatives(mapper, time, theta, thetaH, oneMinusThetaH); } /** {@inheritDoc} */ @@ -202,13 +157,15 @@ public abstract class AbstractFieldStepInterpolator computeInterpolatedStateAndDerivatives(FieldEquationsMapper equationsMapper, - T time, T theta, T oneMinusThetaH) + T time, T theta, + T thetaH, T oneMinusThetaH) throws MaxCountExceededException; } diff --git a/src/main/java/org/apache/commons/math4/ode/sampling/FieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/sampling/FieldStepInterpolator.java index b1e08e3ac..a68d09805 100644 --- a/src/main/java/org/apache/commons/math4/ode/sampling/FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/sampling/FieldStepInterpolator.java @@ -18,7 +18,6 @@ package org.apache.commons.math4.ode.sampling; import org.apache.commons.math4.RealFieldElement; -import org.apache.commons.math4.exception.MaxCountExceededException; import org.apache.commons.math4.ode.FieldODEStateAndDerivative; /** This interface represents an interpolator over the last step @@ -74,14 +73,4 @@ public interface FieldStepInterpolator> { */ boolean isForward(); - /** Copy the instance. - *

The copied instance is guaranteed to be independent from the - * original one. Both can be used with different settings for - * interpolated time without any side effect.

- * @return a deep copy of the instance, which can be used independently. - * @exception MaxCountExceededException if the number of functions evaluations is exceeded - * during step finalization - */ - FieldStepInterpolator copy() throws MaxCountExceededException; - } diff --git a/src/test/java/org/apache/commons/math4/ode/ContinuousOutputFieldModelTest.java b/src/test/java/org/apache/commons/math4/ode/ContinuousOutputFieldModelTest.java new file mode 100644 index 000000000..f7bbe5c9c --- /dev/null +++ b/src/test/java/org/apache/commons/math4/ode/ContinuousOutputFieldModelTest.java @@ -0,0 +1,219 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math4.ode; + +import java.util.Random; + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.ode.nonstiff.DormandPrince54FieldIntegrator; +import org.apache.commons.math4.ode.nonstiff.DormandPrince853FieldIntegrator; +import org.apache.commons.math4.ode.sampling.DummyFieldStepInterpolator; +import org.apache.commons.math4.ode.sampling.FieldStepInterpolator; +import org.apache.commons.math4.util.Decimal64Field; +import org.apache.commons.math4.util.FastMath; +import org.apache.commons.math4.util.MathArrays; +import org.apache.commons.math4.util.MathUtils; +import org.junit.Assert; +import org.junit.Test; + +public class ContinuousOutputFieldModelTest { + + @Test + public void testBoundaries() { + doTestBoundaries(Decimal64Field.getInstance()); + } + + private > void doTestBoundaries(final Field field) { + TestFieldProblem3 pb = new TestFieldProblem3(field, field.getZero().add(0.9)); + double minStep = 0; + double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(); + FieldFirstOrderIntegrator integ = new DormandPrince54FieldIntegrator(field, minStep, maxStep, 1.0e-8, 1.0e-8); + integ.addStepHandler(new ContinuousOutputFieldModel()); + integ.integrate(new FieldExpandableODE(pb), pb.getInitialState(), pb.getFinalTime()); + ContinuousOutputFieldModel cm = (ContinuousOutputFieldModel) integ.getStepHandlers().iterator().next(); + cm.getInterpolatedState(pb.getInitialState().getTime().multiply(2).subtract(pb.getFinalTime())); + cm.getInterpolatedState(pb.getFinalTime().multiply(2).subtract(pb.getInitialState().getTime())); + cm.getInterpolatedState(pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5)); + } + + @Test + public void testRandomAccess() { + doTestRandomAccess(Decimal64Field.getInstance()); + } + + private > void doTestRandomAccess(final Field field) { + + TestFieldProblem3 pb = new TestFieldProblem3(field, field.getZero().add(0.9)); + double minStep = 0; + double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(); + FieldFirstOrderIntegrator integ = new DormandPrince54FieldIntegrator(field, minStep, maxStep, 1.0e-8, 1.0e-8); + ContinuousOutputFieldModel cm = new ContinuousOutputFieldModel(); + integ.addStepHandler(cm); + integ.integrate(new FieldExpandableODE(pb), pb.getInitialState(), pb.getFinalTime()); + + Random random = new Random(347588535632l); + T maxError = field.getZero(); + T maxErrorDot = field.getZero(); + for (int i = 0; i < 1000; ++i) { + double r = random.nextDouble(); + T time = pb.getInitialState().getTime().multiply(r).add(pb.getFinalTime().multiply(1.0 - r)); + FieldODEStateAndDerivative interpolated = cm.getInterpolatedState(time); + T[] theoreticalY = pb.computeTheoreticalState(time); + T[] theoreticalYDot = pb.doComputeDerivatives(time, theoreticalY); + T dx = interpolated.getState()[0].subtract(theoreticalY[0]); + T dy = interpolated.getState()[1].subtract(theoreticalY[1]); + T error = dx.multiply(dx).add(dy.multiply(dy)); + maxError = MathUtils.max(maxError, error); + T dxDot = interpolated.getDerivative()[0].subtract(theoreticalYDot[0]); + T dyDot = interpolated.getDerivative()[1].subtract(theoreticalYDot[1]); + T errorDot = dxDot.multiply(dxDot).add(dyDot.multiply(dyDot)); + maxErrorDot = MathUtils.max(maxErrorDot, errorDot); + } + + Assert.assertEquals(0.0, maxError.getReal(), 1.0e-9); + Assert.assertEquals(0.0, maxErrorDot.getReal(), 4.0e-7); + + } + + @Test + public void testModelsMerging() { + doTestModelsMerging(Decimal64Field.getInstance()); + } + + private > void doTestModelsMerging(final Field field) { + + // theoretical solution: y[0] = cos(t), y[1] = sin(t) + FieldFirstOrderDifferentialEquations problem = + new FieldFirstOrderDifferentialEquations() { + public T[] computeDerivatives(T t, T[] y) { + T[] yDot = MathArrays.buildArray(field, 2); + yDot[0] = y[1].negate(); + yDot[1] = y[0]; + return yDot; + } + public int getDimension() { + return 2; + } + public void init(T t0, T[] y0, T finalTime) { + } + }; + + // integrate backward from π to 0; + ContinuousOutputFieldModel cm1 = new ContinuousOutputFieldModel(); + FieldFirstOrderIntegrator integ1 = + new DormandPrince853FieldIntegrator(field, 0, 1.0, 1.0e-8, 1.0e-8); + integ1.addStepHandler(cm1); + T t0 = field.getZero().add(FastMath.PI); + T[] y0 = MathArrays.buildArray(field, 2); + y0[0] = field.getOne().negate(); + y0[1] = field.getZero(); + integ1.integrate(new FieldExpandableODE(problem), + new FieldODEState(t0, y0), + field.getZero()); + + // integrate backward from 2π to π + ContinuousOutputFieldModel cm2 = new ContinuousOutputFieldModel(); + FieldFirstOrderIntegrator integ2 = + new DormandPrince853FieldIntegrator(field, 0, 0.1, 1.0e-12, 1.0e-12); + integ2.addStepHandler(cm2); + t0 = field.getZero().add(2.0 * FastMath.PI); + y0[0] = field.getOne(); + y0[1] = field.getZero(); + integ2.integrate(new FieldExpandableODE(problem), + new FieldODEState(t0, y0), + field.getZero().add(FastMath.PI)); + + // merge the two half circles + ContinuousOutputFieldModel cm = new ContinuousOutputFieldModel(); + cm.append(cm2); + cm.append(new ContinuousOutputFieldModel()); + cm.append(cm1); + + // check circle + Assert.assertEquals(2.0 * FastMath.PI, cm.getInitialTime().getReal(), 1.0e-12); + Assert.assertEquals(0, cm.getFinalTime().getReal(), 1.0e-12); + for (double t = 0; t < 2.0 * FastMath.PI; t += 0.1) { + FieldODEStateAndDerivative interpolated = cm.getInterpolatedState(field.getZero().add(t)); + Assert.assertEquals(FastMath.cos(t), interpolated.getState()[0].getReal(), 1.0e-7); + Assert.assertEquals(FastMath.sin(t), interpolated.getState()[1].getReal(), 1.0e-7); + } + + } + + @Test + public void testErrorConditions() { + doTestErrorConditions(Decimal64Field.getInstance()); + } + + private > void doTestErrorConditions(final Field field) { + ContinuousOutputFieldModel cm = new ContinuousOutputFieldModel(); + cm.handleStep(buildInterpolator(field, 0, 1, new double[] { 0.0, 1.0, -2.0 }), true); + + // dimension mismatch + Assert.assertTrue(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0 })); + + // hole between time ranges + Assert.assertTrue(checkAppendError(field, cm, 10.0, 20.0, new double[] { 0.0, 1.0, -2.0 })); + + // propagation direction mismatch + Assert.assertTrue(checkAppendError(field, cm, 1.0, 0.0, new double[] { 0.0, 1.0, -2.0 })); + + // no errors + Assert.assertFalse(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0, -2.0 })); + + } + + private > boolean checkAppendError(Field field, ContinuousOutputFieldModel cm, + double t0, double t1, double[] y) { + try { + ContinuousOutputFieldModel otherCm = new ContinuousOutputFieldModel(); + otherCm.handleStep(buildInterpolator(field, t0, t1, y), true); + cm.append(otherCm); + } catch(IllegalArgumentException iae) { + return true; // there was an allowable error + } + return false; // no allowable error + } + + private > FieldStepInterpolator buildInterpolator(Field field, + double t0, double t1, double[] y) { + T[] fieldY = MathArrays.buildArray(field, y.length); + for (int i = 0; i < y.length; ++i) { + fieldY[i] = field.getZero().add(y[i]); + } + final FieldODEStateAndDerivative s0 = new FieldODEStateAndDerivative(field.getZero().add(t0), fieldY, fieldY); + final FieldODEStateAndDerivative s1 = new FieldODEStateAndDerivative(field.getZero().add(t1), fieldY, fieldY); + final FieldEquationsMapper mapper = new FieldExpandableODE(new FieldFirstOrderDifferentialEquations() { + public int getDimension() { + return s0.getStateDimension(); + } + public void init(T t0, T[] y0, T finalTime) { + } + public T[] computeDerivatives(T t, T[] y) { + return y; + } + }).getMapper(); + return new DummyFieldStepInterpolator(t1 >= t0, s0, s1, s0, s1, mapper); + } + + public void checkValue(double value, double reference) { + Assert.assertTrue(FastMath.abs(value - reference) < 1.0e-10); + } + +} diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldStepInterpolatorTest.java index 4716b81e7..b670cda75 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldStepInterpolatorTest.java @@ -18,8 +18,6 @@ package org.apache.commons.math4.ode.nonstiff; -import java.lang.reflect.InvocationTargetException; - import org.apache.commons.math4.Field; import org.apache.commons.math4.RealFieldElement; import org.apache.commons.math4.ode.AbstractIntegrator; @@ -38,7 +36,15 @@ import org.junit.Test; public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { protected abstract > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper); + + protected abstract > FieldButcherArrayProvider + createButcherArrayProvider(final Field field); @Test public abstract void interpolationAtBounds(); @@ -140,10 +146,8 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { final double t0, final double[] y0, final double t1) { - RungeKuttaFieldStepInterpolator interpolator = createInterpolator(field, t1 > t0, - new FieldExpandableODE(eqn).getMapper()); // get the Butcher arrays from the field integrator - FieldButcherArrayProvider provider = createButcherArrayProvider(field, interpolator); + FieldButcherArrayProvider provider = createButcherArrayProvider(field); T[][] a = provider.getA(); T[] b = provider.getB(); T[] c = provider.getC(); @@ -156,8 +160,7 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { fieldY[i] = field.getZero().add(y0[i]); } fieldYDotK[0] = eqn.computeDerivatives(t, fieldY); - interpolator.storeState(new FieldODEStateAndDerivative(t, fieldY, fieldYDotK[0])); - interpolator.shift(); + FieldODEStateAndDerivative s0 = new FieldODEStateAndDerivative(t, fieldY, fieldYDotK[0]); // perform one integration step, in order to get consistent derivatives T h = field.getZero().add(t1 - t0); @@ -170,20 +173,20 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { } fieldYDotK[k + 1] = eqn.computeDerivatives(h.multiply(c[k]).add(t0), fieldY); } - interpolator.setSlopes(fieldYDotK); // store state at step end + t = field.getZero().add(t1); for (int i = 0; i < y0.length; ++i) { fieldY[i] = field.getZero().add(y0[i]); for (int s = 0; s < b.length; ++s) { fieldY[i] = fieldY[i].add(h.multiply(b[s].multiply(fieldYDotK[s][i]))); } } - interpolator.storeState(new FieldODEStateAndDerivative(field.getZero().add(t1), - fieldY, - eqn.computeDerivatives(field.getZero().add(t1), fieldY))); + FieldODEStateAndDerivative s1 = new FieldODEStateAndDerivative(t, fieldY, + eqn.computeDerivatives(t, fieldY)); - return interpolator; + return createInterpolator(field, t1 > t0, fieldYDotK, s0, s1, s0, s1, + new FieldExpandableODE(eqn).getMapper()); } @@ -236,10 +239,10 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { } @Override public void computeDerivatives(final double t, final double[] y, final double[] yDot) { - T fieldT = fieldInterpolator.getField().getZero().add(t); - T[] fieldY = MathArrays.buildArray(fieldInterpolator.getField(), y.length); + T fieldT = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(t); + T[] fieldY = MathArrays.buildArray(fieldInterpolator.getCurrentState().getTime().getField(), y.length); for (int i = 0; i < y.length; ++i) { - fieldY[i] = fieldInterpolator.getField().getZero().add(y[i]); + fieldY[i] = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(y[i]); } T[] fieldYDot = eqn.computeDerivatives(fieldT, fieldY); for (int i = 0; i < yDot.length; ++i) { @@ -281,42 +284,6 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest { } - private > FieldButcherArrayProvider - createButcherArrayProvider(final Field field, final RungeKuttaFieldStepInterpolator provider) { - FieldButcherArrayProvider integrator = null; - try { - String interpolatorName = provider.getClass().getName(); - String integratorName = interpolatorName.replaceAll("StepInterpolator", "Integrator"); - @SuppressWarnings("unchecked") - Class> clz = (Class>) Class.forName(integratorName); - try { - integrator = clz.getConstructor(Field.class, RealFieldElement.class). - newInstance(field, field.getOne()); - } catch (NoSuchMethodException nsme) { - try { - integrator = clz.getConstructor(Field.class, - Double.TYPE, Double.TYPE, - Double.TYPE, Double.TYPE). - newInstance(field, 0.001, 1.0, 1.0, 1.0); - } catch (NoSuchMethodException e) { - Assert.fail(e.getLocalizedMessage()); - } - } - - } catch (InvocationTargetException ite) { - Assert.fail(ite.getLocalizedMessage()); - } catch (IllegalAccessException iae) { - Assert.fail(iae.getLocalizedMessage()); - } catch (InstantiationException ie) { - Assert.fail(ie.getLocalizedMessage()); - } catch (ClassNotFoundException cnfe) { - Assert.fail(cnfe.getLocalizedMessage()); - } - - return integrator; - - } - private static class SinCos> implements FieldFirstOrderDifferentialEquations { private final Field field; protected SinCos(final Field field) { diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungKuttaFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungKuttaFieldStepInterpolatorTest.java index d5592006e..527fe25f2 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungKuttaFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungKuttaFieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class ClassicalRungKuttaFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new ClassicalRungeKuttaFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new ClassicalRungeKuttaFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new ClassicalRungeKuttaFieldIntegrator(field, field.getOne()); } @Test @@ -43,7 +57,7 @@ public class ClassicalRungKuttaFieldStepInterpolatorTest extends AbstractRungeKu @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 1.2e-16, 2.3e-16, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 1.2e-16, 3.4e-16, 2.1e-17); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolatorTest.java index 2aeac5490..ee17388db 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class DormandPrince54FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new DormandPrince54FieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new DormandPrince54FieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new DormandPrince54FieldIntegrator(field, 0, 1, 1, 1); } @Test @@ -43,7 +57,7 @@ public class DormandPrince54FieldStepInterpolatorTest extends AbstractRungeKutta @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 2.3e-16, 4.5e-16, 5.6e-17); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 2.3e-16, 5.6e-16, 5.6e-17); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java index bbef7b466..a104867f2 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class DormandPrince853FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new DormandPrince853FieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new DormandPrince853FieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new DormandPrince853FieldIntegrator(field, 0, 1, 1, 1); } @Test diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolatorTest.java index b7fa85eda..0458f5c3b 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolatorTest.java @@ -21,14 +21,27 @@ 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.Decimal64Field; import org.junit.Test; public class EulerFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new EulerFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, FieldEquationsMapper mapper) { + return new EulerFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new EulerFieldIntegrator(field, field.getOne()); } @Test @@ -43,7 +56,7 @@ public class EulerFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepI @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 7.0e-18, 1.0e-50, 1.0e-50, 1.0e-50); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolatorTest.java index 7dcc87451..54f0098e4 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class GillFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new GillFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new GillFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new GillFieldIntegrator(field, field.getOne()); } @Test @@ -43,7 +57,7 @@ public class GillFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepIn @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 3.4e-16, 2.1e-17); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolatorTest.java index 0d9180b93..4bbe6f025 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class HighamHall54FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new HighamHall54FieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new HighamHall54FieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new HighamHall54FieldIntegrator(field, 0, 1, 1, 1); } @Test @@ -43,7 +57,7 @@ public class HighamHall54FieldStepInterpolatorTest extends AbstractRungeKuttaFie @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.3e-16, 1.0e-50, 3.6e-15, 2.1e-16); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolatorTest.java index 4add8ca0d..39ebdb9f9 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class LutherFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new LutherFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new LutherFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new LutherFieldIntegrator(field, field.getOne()); } @Test diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolatorTest.java index 20dee4780..6779c4ade 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class MidpointFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new MidpointFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new MidpointFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new MidpointFieldIntegrator(field, field.getOne()); } @Test @@ -43,7 +57,7 @@ public class MidpointFieldStepInterpolatorTest extends AbstractRungeKuttaFieldSt @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 7.0e-18); } } diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java index 6b956fae2..6891e02be 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java @@ -21,14 +21,28 @@ 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.Decimal64Field; import org.junit.Test; public class ThreeEighthesFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest { protected > RungeKuttaFieldStepInterpolator - createInterpolator(Field field, boolean forward, FieldEquationsMapper mapper) { - return new ThreeEighthesFieldStepInterpolator(field, forward, mapper); + createInterpolator(Field field, boolean forward, T[][] yDotK, + FieldODEStateAndDerivative globalPreviousState, + FieldODEStateAndDerivative globalCurrentState, + FieldODEStateAndDerivative softPreviousState, + FieldODEStateAndDerivative softCurrentState, + FieldEquationsMapper mapper) { + return new ThreeEighthesFieldStepInterpolator(field, forward, yDotK, + globalPreviousState, globalCurrentState, + softPreviousState, softCurrentState, + mapper); + } + + protected > FieldButcherArrayProvider + createButcherArrayProvider(final Field field) { + return new ThreeEighthesFieldIntegrator(field, field.getOne()); } @Test @@ -43,7 +57,7 @@ public class ThreeEighthesFieldStepInterpolatorTest extends AbstractRungeKuttaFi @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.2e-16, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.2e-16, 3.4e-16, 1.4e-17); } } diff --git a/src/test/java/org/apache/commons/math4/ode/sampling/DummyFieldStepInterpolator.java b/src/test/java/org/apache/commons/math4/ode/sampling/DummyFieldStepInterpolator.java new file mode 100644 index 000000000..f5467459d --- /dev/null +++ b/src/test/java/org/apache/commons/math4/ode/sampling/DummyFieldStepInterpolator.java @@ -0,0 +1,55 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math4.ode.sampling; + +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; + +public class DummyFieldStepInterpolator> + extends AbstractFieldStepInterpolator { + + public DummyFieldStepInterpolator(final boolean forward, + final FieldODEStateAndDerivative globalPreviousState, + final FieldODEStateAndDerivative globalCurrentState, + final FieldODEStateAndDerivative softPreviousState, + final FieldODEStateAndDerivative softCurrentState, + final FieldEquationsMapper mapper) { + super(forward, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper); + } + + @Override + protected AbstractFieldStepInterpolator create(final boolean newForward, + final FieldODEStateAndDerivative newGlobalPreviousState, + final FieldODEStateAndDerivative newGlobalCurrentState, + final FieldODEStateAndDerivative newSoftPreviousState, + final FieldODEStateAndDerivative newSoftCurrentState, + final FieldEquationsMapper newMapper) { + return new DummyFieldStepInterpolator(newForward, + newGlobalPreviousState, newGlobalCurrentState, + newSoftPreviousState, newSoftCurrentState, + newMapper); + } + + @Override + protected FieldODEStateAndDerivative computeInterpolatedStateAndDerivatives(FieldEquationsMapper equationsMapper, + T time, T theta, T thetaH, T oneMinusThetaH) { + return getGlobalCurrentState(); + } + +}