Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

One step towards immutable step interpolators. 

The interpolators do not expect anymore the y and yDot arrays to
be shared with integrator and be updated by it.
This commit is contained in:
Luc Maisonobe 2016-01-06 12:40:23 +01:00
parent aab178594f
commit 6d1fb4dc3e
27 changed files with 583 additions and 1032 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -101,10 +100,8 @@ public class ClassicalRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ClassicalRungeKuttaFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(this, forward, mapper);
}
}

60
src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
View File

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

11
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.MathUtils;
@ -93,7 +92,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
super(field, METHOD_NAME, true,
super(field, METHOD_NAME, 6,
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
e1 = fraction( 71, 57600);
e3 = fraction( -71, 16695);
@ -119,7 +118,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
super(field, METHOD_NAME, true,
super(field, METHOD_NAME, 6,
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
e1 = fraction( 71, 57600);
e3 = fraction( -71, 16695);
@ -190,10 +189,8 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected DormandPrince54FieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(this, forward, mapper);
}
/** {@inheritDoc} */

274
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
View File

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

74
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java
View File

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

616
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
View File

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

31
src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
View File

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

7
src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -86,10 +85,8 @@ public class EulerFieldIntegrator<T extends RealFieldElement<T>> extends RungeKu
/** {@inheritDoc} */
@Override
protected EulerFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(this, forward, mapper);
}
}

28
src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java
View File

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

7
src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -111,10 +110,8 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected GillFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(this, forward, mapper);
}
}

60
src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java
View File

@ -22,7 +22,6 @@ import org.apache.commons.math4.ode.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathArrays;
/**
* This class implements a step interpolator for the Gill fourth
@ -68,17 +67,13 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
GillFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(rkIntegrator, y, yDotArray, forward, mapper);
super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@ -98,6 +93,7 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@ -111,48 +107,30 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot2 = cDot23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeffDot3 = cDot23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeffDot4 = theta.multiply(twoTheta.subtract(1));
final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if ((previousState != null) && (theta.getReal() <= 0.5)) {
final T s = theta.multiply(h).divide(6.0);
final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
interpolatedState[i] = previousState[i].
add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)).
add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
}
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(6.0);
final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
} else {
final T s = oneMinusThetaH.divide(6.0);
final T c23 = s .multiply(twoTheta.add(2).subtract(fourTheta2));
final T s = oneMinusThetaH.divide(-6.0);
final T c23 = s.multiply(twoTheta.add(2).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(5)).add(1));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff4 = s.multiply(fourTheta2.add(theta).add(1));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
interpolatedState[i] = currentState[i].
subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)).
subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
}
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}

11
src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.MathUtils;
@ -64,7 +63,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
super(field, METHOD_NAME, false,
super(field, METHOD_NAME, -1,
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
e = MathArrays.buildArray(field, 7);
e[0] = fraction(-1, 20);
@ -92,7 +91,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
super(field, METHOD_NAME, false,
super(field, METHOD_NAME, -1,
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
e = MathArrays.buildArray(field, 7);
e[0] = fraction(-1, 20);
@ -165,10 +164,8 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected HighamHall54FieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(this, forward, mapper);
}
/** {@inheritDoc} */

61
src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
* This class represents an interpolator over the last step during an
@ -38,17 +37,13 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
HighamHall54FieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(rkIntegrator, y, yDotArray, forward, mapper);
super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@ -68,72 +63,44 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T bDot0 = theta.multiply(theta.multiply(theta.multiply( -10.0 ).add( 16.0 )).add(-15.0 / 2.0)).add(1);
final T bDot1 = time.getField().getZero();
final T bDot2 = theta.multiply(theta.multiply(theta.multiply( 135.0 / 2.0).add(-729.0 / 8.0)).add(459.0 / 16.0));
final T bDot3 = theta.multiply(theta.multiply(theta.multiply(-120.0 ).add( 152.0 )).add(-44.0 ));
final T bDot4 = theta.multiply(theta.multiply(theta.multiply( 125.0 / 2.0).add(-625.0 / 8.0)).add(375.0 / 16.0));
final T bDot5 = theta.multiply( 5.0 / 8.0).multiply(theta.multiply(2).subtract(1));
final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if ((previousState != null) && (theta.getReal() <= 0.5)) {
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T hTheta = h.multiply(theta);
final T b0 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply( -5.0 / 2.0).add( 16.0 / 3.0)).add(-15.0 / 4.0)).add(1));
final T b1 = time.getField().getZero();
final T b2 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0)));
final T b3 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 )));
final T b4 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0)));
final T b5 = hTheta.multiply(theta.multiply(theta.multiply( 5.0 / 12.0)).add( -5.0 / 16.0));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot0 = yDotK[0][i];
final T yDot2 = yDotK[2][i];
final T yDot3 = yDotK[3][i];
final T yDot4 = yDotK[4][i];
final T yDot5 = yDotK[5][i];
interpolatedState[i] = previousState[i].
add(b0.multiply(yDot0)).
add(b2.multiply(yDot2)).
add(b3.multiply(yDot3)).
add(b4.multiply(yDot4)).
add(b5.multiply(yDot5));
interpolatedDerivatives[i] = bDot0.multiply(yDot0).
add(bDot2.multiply(yDot2)).
add(bDot3.multiply(yDot3)).
add(bDot4.multiply(yDot4)).
add(bDot5.multiply(yDot5));
}
interpolatedState = previousStateLinearCombination(b0, b1, b2, b3, b4, b5);
interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5);
} else {
final T theta2 = theta.multiply(theta);
final T b0 = h.multiply( theta.multiply(theta.multiply(theta.multiply(theta.multiply(-5.0 / 2.0).add( 16.0 / 3.0)).add( -15.0 / 4.0)).add( 1.0 )).add( -1.0 / 12.0));
final T b1 = time.getField().getZero();
final T b2 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 135.0 / 8.0 ).add(-243.0 / 8.0)).add(459.0 / 32.0)).add( -27.0 / 32.0));
final T b3 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( -30.0 ).add( 152.0 / 3.0)).add(-22.0 )).add( 4.0 / 3.0));
final T b4 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 125.0 / 8.0 ).add(-625.0 / 24.0)).add(375.0 / 32.0)).add(-125.0 / 96.0));
final T b5 = h.multiply(theta2.multiply(theta.multiply( 5.0 / 12.0 ).add(-5.0 / 16.0)).add( -5.0 / 48.0));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot0 = yDotK[0][i];
final T yDot2 = yDotK[2][i];
final T yDot3 = yDotK[3][i];
final T yDot4 = yDotK[4][i];
final T yDot5 = yDotK[5][i];
interpolatedState[i] = currentState[i].
add(b0.multiply(yDot0)).
add(b2.multiply(yDot2)).
add(b3.multiply(yDot3)).
add(b4.multiply(yDot4)).
add(b5.multiply(yDot5));
interpolatedDerivatives[i] = bDot0.multiply(yDot0).
add(bDot2.multiply(yDot2)).
add(bDot3.multiply(yDot3)).
add(bDot4.multiply(yDot4)).
add(bDot5.multiply(yDot5));
}
interpolatedState = currentStateLinearCombination(b0, b1, b2, b3, b4, b5);
interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5);
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}

7
src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -138,10 +137,8 @@ public class LutherFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected LutherFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(this, forward, mapper);
}
}

104
src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
* This class represents an interpolator over the last step during an
@ -83,17 +82,13 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
LutherFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(rkIntegrator, y, yDotArray, forward, mapper);
super(rkIntegrator, forward, mapper);
final T q = rkIntegrator.getField().getOne().multiply(21).sqrt();
c5a = q.multiply( -49).add( -49);
c5b = q.multiply( 287).add( 392);
@ -143,6 +138,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@ -192,86 +188,42 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
// At the end, we get the b_i as polynomials in theta.
final T coeffDot1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 ).add( -47 )).add( 36 )).add( -54 / 5.0)).add(1);
// not really needed as it is zero: final T coeffDot2 = theta.getField().getZero();
final T coeffDot2 = theta.getField().getZero();
final T coeffDot3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 ).add(-608 / 3.0)).add( 320 / 3.0 )).add(-208 / 15.0));
final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0));
final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150)));
final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150)));
final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3 )).add( -3 )).add( 3 / 5.0);
final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if ((previousState != null) && (theta.getReal() <= 0.5)) {
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1);
// not really needed as it is zero: final T coeff2 = theta.getField().getZero();
final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0));
final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0));
final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300)));
final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300)));
final T coeff7 = theta.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0);
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
// not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
final T yDot5 = yDotK[4][i];
final T yDot6 = yDotK[5][i];
final T yDot7 = yDotK[6][i];
interpolatedState[i] = previousState[i].
add(theta.multiply(h).
multiply( coeff1.multiply(yDot1).
// not really needed as it is zero: add(coeff2.multiply(yDot2)).
add(coeff3.multiply(yDot3)).
add(coeff4.multiply(yDot4)).
add(coeff5.multiply(yDot5)).
add(coeff6.multiply(yDot6)).
add(coeff7.multiply(yDot7))));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).
// not really needed as it is zero: add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).
add(coeffDot4.multiply(yDot4)).
add(coeffDot5.multiply(yDot5)).
add(coeffDot6.multiply(yDot6)).
add(coeffDot7.multiply(yDot7));
}
final T s = theta.multiply(theta.multiply(h));
final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1);
final T coeff2 = s.getField().getZero();
final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0));
final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0));
final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300)));
final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300)));
final T coeff7 = s.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0);
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
} else {
final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0);
// not really needed as it is zero: final T coeff2 = theta.getField().getZero();
final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0);
final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0)));
final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
final T coeff7 = theta.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0);
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
// not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
final T yDot5 = yDotK[4][i];
final T yDot6 = yDotK[5][i];
final T yDot7 = yDotK[6][i];
interpolatedState[i] = currentState[i].
add(oneMinusThetaH.
multiply( coeff1.multiply(yDot1).
// not really needed as it is zero: add(coeff2.multiply(yDot2)).
add(coeff3.multiply(yDot3)).
add(coeff4.multiply(yDot4)).
add(coeff5.multiply(yDot5)).
add(coeff6.multiply(yDot6)).
add(coeff7.multiply(yDot7))));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).
// not really needed as it is zero: add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).
add(coeffDot4.multiply(yDot4)).
add(coeffDot5.multiply(yDot5)).
add(coeffDot6.multiply(yDot6)).
add(coeffDot7.multiply(yDot7));
}
final T s = oneMinusThetaH.multiply(theta);
final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0);
final T coeff2 = s.getField().getZero();
final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0);
final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0)));
final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
final T coeff7 = s.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0);
interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}

7
src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -86,10 +85,8 @@ public class MidpointFieldIntegrator<T extends RealFieldElement<T>> extends Rung
/** {@inheritDoc} */
@Override
protected MidpointFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(this, forward, mapper);
}
}

40
src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
* This class implements a step interpolator for second order
@ -54,17 +53,13 @@ class MidpointFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
MidpointFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(rkIntegrator, y, yDotArray, forward, mapper);
super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@ -84,37 +79,30 @@ class MidpointFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T coeffDot2 = theta.multiply(2);
final T coeffDot1 = time.getField().getOne().subtract(coeffDot2);
final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
final T coeffDot2 = theta.multiply(2);
final T coeffDot1 = time.getField().getOne().subtract(coeffDot2);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if ((previousState != null) && (theta.getReal() <= 0.5)) {
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T coeff1 = theta.multiply(oneMinusThetaH);
final T coeff2 = theta.multiply(theta).multiply(h);
for (int i = 0; i < previousState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
interpolatedState[i] = previousState[i].add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2));
}
interpolatedState = previousStateLinearCombination(coeff1, coeff2);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2);
} else {
final T coeff1 = oneMinusThetaH.multiply(theta);
final T coeff2 = oneMinusThetaH.multiply(theta.add(1));
for (int i = 0; i < previousState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
interpolatedState[i] = currentState[i].add(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2));
}
final T coeff2 = oneMinusThetaH.multiply(theta.add(1)).negate();
interpolatedState = currentStateLinearCombination(coeff1, coeff2);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2);
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}

13
src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java
View File

@ -112,18 +112,11 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
protected abstract T[] getB();
/** Create an interpolator.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
* @return external weights for the high order method from Butcher array
*/
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(AbstractFieldIntegrator<T> rkIntegrator,
T[] y, T[][] yDotArray,
boolean forward,
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward,
FieldEquationsMapper<T> mapper);
/** {@inheritDoc} */
@ -147,7 +140,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
// set up an interpolator sharing the integrator arrays
final RungeKuttaFieldStepInterpolator<T> interpolator =
createInterpolator(this, y0, yDotK, forward, equations.getMapper());
createInterpolator(forward, equations.getMapper());
interpolator.storeState(stepStart);
// set up integration control objects
@ -172,6 +165,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
interpolator.shift();
// first stage
y = equations.getMapper().mapState(stepStart);
yDotK[0] = equations.getMapper().mapDerivative(stepStart);
// next stages
@ -202,6 +196,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
// discrete events handling
interpolator.setSlopes(yDotK);
interpolator.storeState(stateTmp);
System.arraycopy(yTmp, 0, y, 0, y0.length);
stepStart = acceptStep(interpolator, finalTime);

89
src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.java
View File

@ -36,31 +36,23 @@ import org.apache.commons.math4.util.MathArrays;
abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
extends AbstractFieldStepInterpolator<T> {
/** Previous state. */
protected T[] previousState;
/** Slopes at the intermediate points */
protected T[][] yDotK;
/** Reference to the integrator. */
protected AbstractFieldIntegrator<T> integrator;
/** Slopes at the intermediate points. */
private T[][] yDotK;
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
protected RungeKuttaFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(y, forward, mapper);
this.previousState = null;
this.yDotK = yDotArray;
this.integrator = rkIntegrator;
super(forward, mapper);
this.yDotK = null;
this.integrator = rkIntegrator;
}
/** Copy constructor.
@ -84,18 +76,14 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
super(interpolator);
if (interpolator.currentState != null) {
previousState = interpolator.previousState.clone();
if (yDotK != null) {
yDotK = MathArrays.buildArray(interpolator.integrator.getField(),
interpolator.yDotK.length, interpolator.yDotK[0].length);
for (int k = 0; k < interpolator.yDotK.length; ++k) {
System.arraycopy(interpolator.yDotK[k], 0, yDotK[k], 0, interpolator.yDotK[k].length);
interpolator.yDotK.length, -1);
for (int k = 0; k < yDotK.length; ++k) {
yDotK[k] = interpolator.yDotK[k].clone();
}
} else {
previousState = null;
yDotK = null;
}
@ -105,11 +93,56 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
}
/** {@inheritDoc} */
@Override
public void shift() {
previousState = currentState.clone();
super.shift();
/** Store the slopes at the intermediate points.
* @param slopes slopes at the intermediate points
*/
void setSlopes(final T[][] slopes) {
this.yDotK = slopes.clone();
}
/** Compute a state by linear combination added to previous state.
* @param coefficients coefficients to apply to the method staged derivatives
* @return combined state
*/
@SafeVarargs
protected final T[] previousStateLinearCombination(final T ... coefficients) {
return combine(getPreviousState().getState(),
coefficients);
}
/** Compute a state by linear combination added to current state.
* @param coefficients coefficients to apply to the method staged derivatives
* @return combined state
*/
@SuppressWarnings("unchecked")
protected T[] currentStateLinearCombination(final T ... coefficients) {
return combine(getCurrentState().getState(),
coefficients);
}
/** Compute a state derivative by linear combination.
* @param coefficients coefficients to apply to the method staged derivatives
* @return combined state
*/
@SuppressWarnings("unchecked")
protected T[] derivativeLinearCombination(final T ... coefficients) {
return combine(MathArrays.buildArray(integrator.getField(), yDotK[0].length),
coefficients);
}
/** Linearly combine arrays.
* @param a array to add to
* @param coefficients coefficients to apply to the method staged derivatives
* @return a itself, as a conveniency for fluent API
*/
@SuppressWarnings("unchecked")
private T[] combine(final T[] a, final T ... coefficients) {
for (int i = 0; i < a.length; ++i) {
for (int k = 0; k < coefficients.length; ++k) {
a[i] = a[i].add(coefficients[k].multiply(yDotK[k][i]));
}
}
return a;
}
}

7
src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java
View File

@ -19,7 +19,6 @@ 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.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.util.MathArrays;
@ -100,10 +99,8 @@ public class ThreeEighthesFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ThreeEighthesFieldStepInterpolator<T>
createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
final T[][] yDotArray, final boolean forward,
final FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(this, forward, mapper);
}
}

48
src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.AbstractFieldIntegrator;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
* This class implements a step interpolator for the 3/8 fourth
@ -64,17 +63,13 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
ThreeEighthesFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
final T[] y, final T[][] yDotArray, final boolean forward,
final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(rkIntegrator, y, yDotArray, forward, mapper);
super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@ -94,6 +89,7 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@ -103,51 +99,31 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1);
final T coeffDot2 = coeffDot3.multiply(theta.multiply(-6).add(5));
final T coeffDot4 = coeffDot3.multiply(theta.multiply(2).subtract(1));
final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if ((previousState != null) && (theta.getReal() <= 0.5)) {
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(8);
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(15)).add(8));
final T coeff2 = s.multiply(theta.multiply(5).subtract(fourTheta2)).multiply(3);
final T coeff3 = s.multiply(theta).multiply(3);
final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
interpolatedState[i] = previousState[i].
add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)).
add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
}
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
} else {
final T s = oneMinusThetaH.divide(8);
final T s = oneMinusThetaH.divide(-8);
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T thetaPlus1 = theta.add(1);
final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(7)).add(1));
final T coeff2 = s.multiply(thetaPlus1.subtract(fourTheta2)).multiply(3);
final T coeff3 = s.multiply(thetaPlus1).multiply(3);
final T coeff4 = s.multiply(thetaPlus1.add(fourTheta2));
for (int i = 0; i < interpolatedState.length; ++i) {
final T yDot1 = yDotK[0][i];
final T yDot2 = yDotK[1][i];
final T yDot3 = yDotK[2][i];
final T yDot4 = yDotK[3][i];
interpolatedState[i] = currentState[i].
subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)).
subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4));
interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
}
interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}

21
src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java
View File

@ -43,9 +43,6 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
/** Current time step. */
protected T h;
/** Current state. */
protected T[] currentState;
/** Global previous state. */
private FieldODEStateAndDerivative<T> globalPreviousState;
@ -68,18 +65,16 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
private FieldEquationsMapper<T> mapper;
/** Simple constructor.
* @param y reference to the integrator array holding the state at the end of the step
* @param isForward integration direction indicator
* @param equationsMapper mapper for ODE equations primary and secondary components
*/
protected AbstractFieldStepInterpolator(final T[] y, final boolean isForward,
protected AbstractFieldStepInterpolator(final boolean isForward,
final FieldEquationsMapper<T> equationsMapper) {
globalPreviousState = null;
globalCurrentState = null;
softPreviousState = null;
softCurrentState = null;
h = null;
currentState = y;
finalized = false;
this.forward = isForward;
this.mapper = equationsMapper;
@ -109,17 +104,9 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
softPreviousState = interpolator.softPreviousState;
softCurrentState = interpolator.softCurrentState;
h = interpolator.h;
if (interpolator.currentState == null) {
currentState = null;
mapper = null;
} else {
currentState = interpolator.currentState.clone();
}
finalized = interpolator.finalized;
forward = interpolator.forward;
mapper = interpolator.mapper;
finalized = interpolator.finalized;
forward = interpolator.forward;
mapper = interpolator.mapper;
}

5
src/test/java/org/apache/commons/math4/ode/TestFieldProblem1.java
View File

@ -66,9 +66,10 @@ public class TestFieldProblem1<T extends RealFieldElement<T>>
@Override
public T[] computeTheoreticalState(T t) {
final T[] y0 = getInitialState();
final FieldODEState<T> s0 = getInitialState();
final T[] y0 = s0.getState();
final T[] y = MathArrays.buildArray(getField(), getDimension());
T c = getInitialTime().subtract(t).exp();
T c = s0.getTime().subtract(t).exp();
for (int i = 0; i < getDimension(); ++i) {
y[i] = c.multiply(y0[i]);
}

2
src/test/java/org/apache/commons/math4/ode/TestFieldProblem5.java
View File

@ -35,7 +35,7 @@ public class TestFieldProblem5<T extends RealFieldElement<T>>
*/
public TestFieldProblem5(Field<T> field) {
super(field);
setFinalConditions(getInitialTime().multiply(2).subtract(getFinalTime()));
setFinalConditions(getInitialState().getTime().multiply(2).subtract(getFinalTime()));
}
}

18
src/test/java/org/apache/commons/math4/ode/TestFieldProblemAbstract.java
View File

@ -109,20 +109,12 @@ public abstract class TestFieldProblemAbstract<T extends RealFieldElement<T>>
return n;
}
/**
* Get the initial time.
* @return initial time
/**
* Get the initial state.
* @return initial state
*/
public T getInitialTime() {
return t0;
}
/**
* Get the initial state vector.
* @return initial state vector
*/
public T[] getInitialState() {
return y0;
public FieldODEState<T> getInitialState() {
return new FieldODEState<T>(t0, y0);
}
/**

2
src/test/java/org/apache/commons/math4/ode/TestFieldProblemHandler.java
View File

@ -74,7 +74,7 @@ public class TestFieldProblemHandler<T extends RealFieldElement<T>>
public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) throws MaxCountExceededException {
T start = integrator.getCurrentStepStart().getTime();
if (start.subtract(problem.getInitialTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) {
if (start.subtract(problem.getInitialState().getTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) {
// multistep integrators do not handle the first steps themselves
// so we have to make sure the integrator we look at has really started its work
if (expectedStepStart != null) {

5
src/test/java/org/apache/commons/math4/ode/nonstiff/EulerIntegratorTest.java
View File

@ -59,8 +59,8 @@ public class EulerIntegratorTest {
MaxCountExceededException, NoBracketingException {
for (TestProblemAbstract pb : new TestProblemAbstract[] {
new TestProblem1(), new TestProblem2(), new TestProblem3(),
new TestProblem4(), new TestProblem5(), new TestProblem6()
//new TestProblem1(), new TestProblem2(), new TestProblem3(),
new TestProblem4()//, new TestProblem5(), new TestProblem6()
}) {
double previousValueError = Double.NaN;
@ -95,6 +95,7 @@ public class EulerIntegratorTest {
}
previousTimeError = timeError;
System.exit(1);
}
}