diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java index 6487ccf51..68573fd76 100644 --- a/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java +++ b/src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java @@ -55,7 +55,7 @@ class DormandPrince853FieldStepInterpolator> d = MathArrays.buildArray(getField(), 7, 16); // this row is the same as the b array - d[0][ 0] = fraction(104257, 1929240); + d[0][ 0] = fraction(104257, 1920240); d[0][ 1] = getField().getZero(); d[0][ 2] = getField().getZero(); d[0][ 3] = getField().getZero(); @@ -101,10 +101,10 @@ class DormandPrince853FieldStepInterpolator> 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][12] = d[0][12].multiply(2).subtract(1); // really -1 + d[2][13] = d[0][13].multiply(2); // really 0 + d[2][14] = d[0][14].multiply(2); // really 0 + d[2][15] = d[0][15].multiply(2); // really 0 d[3][ 0] = fraction( -17751989329.0, 2106076560.0); d[3][ 1] = getField().getZero(); @@ -215,7 +215,7 @@ class DormandPrince853FieldStepInterpolator> final T oneMinusThetaH) throws MaxCountExceededException { - final T one = theta.getField().getOne(); + final T one = getField().getOne(); final T eta = one.subtract(theta); final T twoTheta = theta.multiply(2); final T theta2 = theta.multiply(theta); @@ -228,6 +228,7 @@ class DormandPrince853FieldStepInterpolator> final T[] interpolatedState; final T[] interpolatedDerivatives; + if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { final T f0 = theta.multiply(h); final T f1 = f0.multiply(eta); @@ -236,18 +237,58 @@ class DormandPrince853FieldStepInterpolator> 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); + final T[] p = MathArrays.buildArray(getField(), 16); + final T[] q = MathArrays.buildArray(getField(), 16); + for (int i = 0; i < p.length; ++i) { + p[i] = f0.multiply(d[0][i]). + add(f1.multiply(d[1][i])). + add(f2.multiply(d[2][i])). + add(f3.multiply(d[3][i])). + add(f4.multiply(d[4][i])). + add(f5.multiply(d[5][i])). + add(f6.multiply(d[6][i])); + q[i] = d[0][i]. + add(dot1.multiply(d[1][i])). + add(dot2.multiply(d[2][i])). + add(dot3.multiply(d[3][i])). + add(dot4.multiply(d[4][i])). + add(dot5.multiply(d[5][i])). + add(dot6.multiply(d[6][i])); + } + interpolatedState = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7], + p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]); + interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7], + q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]); } else { - final T f0 = oneMinusThetaH; + final T f0 = oneMinusThetaH.negate(); 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); + final T[] p = MathArrays.buildArray(getField(), 16); + final T[] q = MathArrays.buildArray(getField(), 16); + for (int i = 0; i < p.length; ++i) { + p[i] = f0.multiply(d[0][i]). + add(f1.multiply(d[1][i])). + add(f2.multiply(d[2][i])). + add(f3.multiply(d[3][i])). + add(f4.multiply(d[4][i])). + add(f5.multiply(d[5][i])). + add(f6.multiply(d[6][i])); + q[i] = d[0][i]. + add(dot1.multiply(d[1][i])). + add(dot2.multiply(d[2][i])). + add(dot3.multiply(d[3][i])). + add(dot4.multiply(d[4][i])). + add(dot5.multiply(d[5][i])). + add(dot6.multiply(d[6][i])); + } + interpolatedState = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7], + p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]); + interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7], + q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]); } return new FieldODEStateAndDerivative(time, interpolatedState, interpolatedDerivatives); diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java index 357d92a55..3f2b0becb 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java @@ -64,6 +64,20 @@ public abstract class AbstractEmbeddedRungeKuttaFieldIntegratorTest { T[][] fieldA = fieldIntegrator.getA(); T[] fieldB = fieldIntegrator.getB(); T[] fieldC = fieldIntegrator.getC(); + if (fieldIntegrator instanceof DormandPrince853FieldIntegrator) { + // special case for Dormand-Prince 8(5,3), the array in the regular + // integrator is smaller because as of 3.X, the interpolation steps + // are not performed by the integrator itself + T[][] reducedFieldA = MathArrays.buildArray(field, 12, -1); + T[] reducedFieldB = MathArrays.buildArray(field, 13); + T[] reducedFieldC = MathArrays.buildArray(field, 12); + System.arraycopy(fieldA, 0, reducedFieldA, 0, reducedFieldA.length); + System.arraycopy(fieldB, 0, reducedFieldB, 0, reducedFieldB.length); + System.arraycopy(fieldC, 0, reducedFieldC, 0, reducedFieldC.length); + fieldA = reducedFieldA; + fieldB = reducedFieldB; + fieldC = reducedFieldC; + } String fieldName = fieldIntegrator.getClass().getName(); String regularName = fieldName.replaceAll("Field", ""); diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegratorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegratorTest.java index 47a59ba5c..ca777ba5e 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegratorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegratorTest.java @@ -49,12 +49,12 @@ public class DormandPrince853FieldIntegratorTest extends AbstractEmbeddedRungeKu @Test public void testBackward() { - doTestBackward(Decimal64Field.getInstance(), 1.1e-7, 1.1e-7, 1.0e-12, "Dormand-Prince 8 (5, 3)"); + doTestBackward(Decimal64Field.getInstance(), 8.1e-8, 1.1e-7, 1.0e-12, "Dormand-Prince 8 (5, 3)"); } @Test public void testKepler() { - doTestKepler(Decimal64Field.getInstance(), 2.4e-10); + doTestKepler(Decimal64Field.getInstance(), 4.4e-11); } @Override diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java index 1aec28402..bbef7b466 100644 --- a/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java @@ -38,12 +38,12 @@ public class DormandPrince853FieldStepInterpolatorTest extends AbstractRungeKutt @Test public void interpolationInside() { - doInterpolationInside(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50); + doInterpolationInside(Decimal64Field.getInstance(), 3.1e-17, 3.4e-16); } @Test public void nonFieldInterpolatorConsistency() { - doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50, 1.0e-50, 1.0e-50); + doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 3.4e-12, 5.7e-11, 1.9e-10, 3.1e-9); } }