commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject [46/50] [abbrv] [math] Fixed field-based Dormand-Prince 8(5, 3) step interpolator.
Date Wed, 09 Dec 2015 16:17:26 GMT
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/MATH_3_X
Commit: 01026aa09dc1b7394975fd9dca41318054cd4216
Parents: cb51342
Author: Luc Maisonobe <luc@apache.org>
Authored: Wed Dec 9 15:56:36 2015 +0100
Committer: Luc Maisonobe <luc@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);
     }
 
 }


Mime
View raw message