You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2015/12/09 17:01:14 UTC
[17/21] [math] Fixed field-based Dormand-Prince 8(5,
3) step interpolator.
Fixed field-based Dormand-Prince 8(5,3) step interpolator.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/01026aa0
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/01026aa0
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/01026aa0
Branch: refs/heads/field-ode
Commit: 01026aa09dc1b7394975fd9dca41318054cd4216
Parents: cb51342
Author: Luc Maisonobe <lu...@apache.org>
Authored: Wed Dec 9 15:56:36 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Wed Dec 9 15:57:11 2015 +0100
----------------------------------------------------------------------
.../DormandPrince853FieldStepInterpolator.java | 63 ++++++++++++++++----
...ctEmbeddedRungeKuttaFieldIntegratorTest.java | 14 +++++
.../DormandPrince853FieldIntegratorTest.java | 4 +-
...rmandPrince853FieldStepInterpolatorTest.java | 4 +-
4 files changed, 70 insertions(+), 15 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/01026aa0/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
index 211f18c..4890847 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
@@ -55,7 +55,7 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
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<T extends RealFieldElement<T>>
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<T extends RealFieldElement<T>>
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<T extends RealFieldElement<T>>
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<T extends RealFieldElement<T>>
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<T>(time, interpolatedState, interpolatedDerivatives);
http://git-wip-us.apache.org/repos/asf/commons-math/blob/01026aa0/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
index 06a05c9..e166f49 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/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", "");
http://git-wip-us.apache.org/repos/asf/commons-math/blob/01026aa0/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegratorTest.java
index 61a9fd8..4932580 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/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
http://git-wip-us.apache.org/repos/asf/commons-math/blob/01026aa0/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java
index 8153ed3..2fd7ca5 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/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);
}
}