Converted constants for embedded Runge-Kutta integrators. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/b051dbda Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/b051dbda Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/b051dbda
Branch: refs/heads/MATH_3_X Commit: b051dbda7cef3c94196f767678de9e3aaeb0d96d Parents: 1aad860 Author: Luc Maisonobe <l...@apache.org> Authored: Sun Nov 29 18:09:21 2015 +0100 Committer: Luc Maisonobe <l...@apache.org> Committed: Sun Nov 29 18:09:21 2015 +0100 ---------------------------------------------------------------------- .../DormandPrince54FieldIntegrator.java | 132 +++++-- .../DormandPrince853FieldIntegrator.java | 358 ++++++++++++------- .../nonstiff/HighamHall54FieldIntegrator.java | 122 +++++-- 3 files changed, 418 insertions(+), 194 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/b051dbda/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java index 7af0375..7475c66 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java @@ -19,6 +19,9 @@ package org.apache.commons.math3.ode.nonstiff; import org.apache.commons.math3.Field; import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.util.MathArrays; import org.apache.commons.math3.util.MathUtils; @@ -54,45 +57,25 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> /** Integrator method name. */ private static final String METHOD_NAME = "Dormand-Prince 5(4)"; - /** Time steps Butcher array. */ - private static final double[] STATIC_C = { - 1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0 - }; - - /** Internal weights Butcher array. */ - private static final double[][] STATIC_A = { - {1.0/5.0}, - {3.0/40.0, 9.0/40.0}, - {44.0/45.0, -56.0/15.0, 32.0/9.0}, - {19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0}, - {9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0}, - {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0} - }; - - /** Propagation weights Butcher array. */ - private static final double[] STATIC_B = { - 35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0 - }; - /** Error array, element 1. */ - private static final double E1 = 71.0 / 57600.0; + private final T e1; // element 2 is zero, so it is neither stored nor used /** Error array, element 3. */ - private static final double E3 = -71.0 / 16695.0; + private final T e3; /** Error array, element 4. */ - private static final double E4 = 71.0 / 1920.0; + private final T e4; /** Error array, element 5. */ - private static final double E5 = -17253.0 / 339200.0; + private final T e5; /** Error array, element 6. */ - private static final double E6 = 22.0 / 525.0; + private final T e6; /** Error array, element 7. */ - private static final double E7 = -1.0 / 40.0; + private final T e7; /** Simple constructor. * Build a fifth order Dormand-Prince integrator with the given step bounds @@ -110,9 +93,14 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { - super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, - new DormandPrince54FieldStepInterpolator<T>(), + super(field, METHOD_NAME, true, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + e1 = fraction( 71, 57600); + e3 = fraction( -71, 16695); + e4 = fraction( 71, 1920); + e5 = fraction(-17253, 339200); + e6 = fraction( 22, 525); + e7 = fraction( -1, 40); } /** Simple constructor. @@ -131,9 +119,81 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { - super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, - new DormandPrince54FieldStepInterpolator<T>(), + super(field, METHOD_NAME, true, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + e1 = fraction( 71, 57600); + e3 = fraction( -71, 16695); + e4 = fraction( 71, 1920); + e5 = fraction(-17253, 339200); + e6 = fraction( 22, 525); + e7 = fraction( -1, 40); + } + + /** {@inheritDoc} */ + @Override + protected T[] getC() { + final T[] c = MathArrays.buildArray(getField(), 6); + c[0] = fraction(1, 5); + c[1] = fraction(3, 10); + c[2] = fraction(5, 5); + c[3] = fraction(8, 9); + c[4] = getField().getOne(); + c[5] = getField().getOne(); + return c; + } + + /** {@inheritDoc} */ + @Override + protected T[][] getA() { + final T[][] a = MathArrays.buildArray(getField(), 6, -1); + for (int i = 0; i < a.length; ++i) { + a[i] = MathArrays.buildArray(getField(), i + 1); + } + a[0][0] = fraction( 1, 5); + a[1][0] = fraction( 3, 4); + a[1][1] = fraction( 9, 40); + a[2][0] = fraction( 44, 45); + a[2][1] = fraction( -56, 15); + a[2][2] = fraction( 32, 9); + a[3][0] = fraction( 19372, 6561); + a[3][1] = fraction(-25360, 2187); + a[3][2] = fraction( 64448, 6561); + a[3][3] = fraction( -212, 729); + a[4][0] = fraction( 9017, 3168); + a[4][1] = fraction( -355, 33); + a[4][2] = fraction( 46732, 5247); + a[4][3] = fraction( 49, 176); + a[4][4] = fraction( -5103, 18656); + a[5][0] = fraction( 35, 384); + a[5][1] = getField().getZero(); + a[5][2] = fraction( 500, 1113); + a[5][3] = fraction( 125, 192); + a[5][4] = fraction( -2187, 6784); + a[5][5] = fraction( 11, 84); + return a; + } + + /** {@inheritDoc} */ + @Override + protected T[] getB() { + final T[] b = MathArrays.buildArray(getField(), 7); + b[0] = fraction( 35, 384); + b[1] = getField().getZero(); + b[2] = fraction( 500, 1113); + b[3] = fraction( 125, 192); + b[4] = fraction(-2187, 6784); + b[5] = fraction( 11, 84); + b[6] = getField().getZero(); + return b; + } + + /** {@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); } /** {@inheritDoc} */ @@ -149,12 +209,12 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>> T error = getField().getZero(); for (int j = 0; j < mainSetDimension; ++j) { - final T errSum = yDotK[0][j].multiply(E1). - add(yDotK[2][j].multiply(E3)). - add(yDotK[3][j].multiply(E4)). - add(yDotK[4][j].multiply(E5)). - add(yDotK[5][j].multiply(E6)). - add(yDotK[6][j].multiply(E7)); + final T errSum = yDotK[0][j].multiply(e1). + add(yDotK[2][j].multiply(e3)). + add(yDotK[3][j].multiply(e4)). + add(yDotK[4][j].multiply(e5)). + add(yDotK[5][j].multiply(e6)). + add(yDotK[6][j].multiply(e7)); final T yScale = MathUtils.max(y0[j].abs(), y1[j].abs()); final T tol = (vecAbsoluteTolerance == null) ? http://git-wip-us.apache.org/repos/asf/commons-math/blob/b051dbda/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java index c4a574c..6f6b436 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java @@ -19,7 +19,9 @@ package org.apache.commons.math3.ode.nonstiff; import org.apache.commons.math3.Field; import org.apache.commons.math3.RealFieldElement; -import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.util.MathArrays; import org.apache.commons.math3.util.MathUtils; @@ -63,149 +65,58 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> /** Integrator method name. */ private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)"; - /** Time steps Butcher array. */ - private static final double[] STATIC_C = { - (12.0 - 2.0 * FastMath.sqrt(6.0)) / 135.0, (6.0 - FastMath.sqrt(6.0)) / 45.0, (6.0 - FastMath.sqrt(6.0)) / 30.0, - (6.0 + FastMath.sqrt(6.0)) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0, - 6.0/7.0, 1.0, 1.0 - }; - - /** Internal weights Butcher array. */ - private static final double[][] STATIC_A = { - - // k2 - {(12.0 - 2.0 * FastMath.sqrt(6.0)) / 135.0}, - - // k3 - {(6.0 - FastMath.sqrt(6.0)) / 180.0, (6.0 - FastMath.sqrt(6.0)) / 60.0}, - - // k4 - {(6.0 - FastMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - FastMath.sqrt(6.0)) / 40.0}, - - // k5 - {(462.0 + 107.0 * FastMath.sqrt(6.0)) / 3000.0, 0.0, - (-402.0 - 197.0 * FastMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * FastMath.sqrt(6.0)) / 375.0}, - - // k6 - {1.0 / 27.0, 0.0, 0.0, (16.0 + FastMath.sqrt(6.0)) / 108.0, (16.0 - FastMath.sqrt(6.0)) / 108.0}, - - // k7 - {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * FastMath.sqrt(6.0)) / 1024.0, - (118.0 - 23.0 * FastMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0}, - - // k8 - {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * FastMath.sqrt(6.0)) / 371293.0, - (51544.0 - 4784.0 * FastMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0}, - - // k9 - {58656157643.0 / 93983540625.0, 0.0, 0.0, - (-1324889724104.0 - 318801444819.0 * FastMath.sqrt(6.0)) / 626556937500.0, - (-1324889724104.0 + 318801444819.0 * FastMath.sqrt(6.0)) / 626556937500.0, - 96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0, - -165125654.0 / 3796875.0}, - - // k10 - {8909899.0 / 18653125.0, 0.0, 0.0, - (-4521408.0 - 1137963.0 * FastMath.sqrt(6.0)) / 2937500.0, - (-4521408.0 + 1137963.0 * FastMath.sqrt(6.0)) / 2937500.0, - 96663078.0 / 4553125.0, 2107245056.0 / 137915625.0, - -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0}, - - // k11 - {-20401265806.0 / 21769653311.0, 0.0, 0.0, - (354216.0 + 94326.0 * FastMath.sqrt(6.0)) / 112847.0, - (354216.0 - 94326.0 * FastMath.sqrt(6.0)) / 112847.0, - -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0, - 14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0, - -1477884375.0 / 485066827.0}, - - // k12 - {39815761.0 / 17514443.0, 0.0, 0.0, - (-3457480.0 - 960905.0 * FastMath.sqrt(6.0)) / 551636.0, - (-3457480.0 + 960905.0 * FastMath.sqrt(6.0)) / 551636.0, - -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0, - -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0, - 226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0}, - - // k13 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 - {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0, - 66578432.0/35198415.0, -1674902723.0/288716400.0, - 54980371265625.0/176692375811392.0, -734375.0/4826304.0, - 171414593.0/851261400.0, 137909.0/3084480.0} - - }; - - /** Propagation weights Butcher array. */ - private static final double[] STATIC_B = { - 104257.0/1920240.0, - 0.0, - 0.0, - 0.0, - 0.0, - 3399327.0/763840.0, - 66578432.0/35198415.0, - -1674902723.0/288716400.0, - 54980371265625.0/176692375811392.0, - -734375.0/4826304.0, - 171414593.0/851261400.0, - 137909.0/3084480.0, - 0.0 - }; - /** First error weights array, element 1. */ - private static final double E1_01 = 116092271.0 / 8848465920.0; + private final T e1_01; // elements 2 to 5 are zero, so they are neither stored nor used /** First error weights array, element 6. */ - private static final double E1_06 = -1871647.0 / 1527680.0; + private final T e1_06; /** First error weights array, element 7. */ - private static final double E1_07 = -69799717.0 / 140793660.0; + private final T e1_07; /** First error weights array, element 8. */ - private static final double E1_08 = 1230164450203.0 / 739113984000.0; + private final T e1_08; /** First error weights array, element 9. */ - private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0; + private final T e1_09; /** First error weights array, element 10. */ - private static final double E1_10 = 464500805.0 / 1389975552.0; + private final T e1_10; /** First error weights array, element 11. */ - private static final double E1_11 = 1606764981773.0 / 19613062656000.0; + private final T e1_11; /** First error weights array, element 12. */ - private static final double E1_12 = -137909.0 / 6168960.0; + private final T e1_12; /** Second error weights array, element 1. */ - private static final double E2_01 = -364463.0 / 1920240.0; + private final T e2_01; // elements 2 to 5 are zero, so they are neither stored nor used /** Second error weights array, element 6. */ - private static final double E2_06 = 3399327.0 / 763840.0; + private final T e2_06; /** Second error weights array, element 7. */ - private static final double E2_07 = 66578432.0 / 35198415.0; + private final T e2_07; /** Second error weights array, element 8. */ - private static final double E2_08 = -1674902723.0 / 288716400.0; + private final T e2_08; /** Second error weights array, element 9. */ - private static final double E2_09 = -74684743568175.0 / 176692375811392.0; + private final T e2_09; /** Second error weights array, element 10. */ - private static final double E2_10 = -734375.0 / 4826304.0; + private final T e2_10; /** Second error weights array, element 11. */ - private static final double E2_11 = 171414593.0 / 851261400.0; + private final T e2_11; /** Second error weights array, element 12. */ - private static final double E2_12 = 69869.0 / 3084480.0; + private final T e2_12; /** Simple constructor. * Build an eighth order Dormand-Prince integrator with the given step bounds @@ -223,9 +134,24 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { - super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, - new DormandPrince853FieldStepInterpolator<T>(), + super(field, METHOD_NAME, true, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + e1_01 = fraction( 116092271.0, 8848465920.0); + e1_06 = fraction( -1871647.0, 1527680.0); + e1_07 = fraction( -69799717.0, 140793660.0); + e1_08 = fraction( 1230164450203.0, 739113984000.0); + e1_09 = fraction(-1980813971228885.0, 5654156025964544.0); + e1_10 = fraction( 464500805.0, 1389975552.0); + e1_11 = fraction( 1606764981773.0, 19613062656000.0); + e1_12 = fraction( -137909.0, 6168960.0); + e2_01 = fraction( -364463.0, 1920240.0); + e2_06 = fraction( 3399327.0, 763840.0); + e2_07 = fraction( 66578432.0, 35198415.0); + e2_08 = fraction( -1674902723.0, 288716400.0); + e2_09 = fraction( -74684743568175.0, 176692375811392.0); + e2_10 = fraction( -734375.0, 4826304.0); + e2_11 = fraction( 171414593.0, 851261400.0); + e2_12 = fraction( 69869.0, 3084480.0); } /** Simple constructor. @@ -244,9 +170,185 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { - super(field, METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, - new DormandPrince853FieldStepInterpolator<T>(), + super(field, METHOD_NAME, true, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + e1_01 = fraction( 116092271.0, 8848465920.0); + e1_06 = fraction( -1871647.0, 1527680.0); + e1_07 = fraction( -69799717.0, 140793660.0); + e1_08 = fraction( 1230164450203.0, 739113984000.0); + e1_09 = fraction(-1980813971228885.0, 5654156025964544.0); + e1_10 = fraction( 464500805.0, 1389975552.0); + e1_11 = fraction( 1606764981773.0, 19613062656000.0); + e1_12 = fraction( -137909.0, 6168960.0); + e2_01 = fraction( -364463.0, 1920240.0); + e2_06 = fraction( 3399327.0, 763840.0); + e2_07 = fraction( 66578432.0, 35198415.0); + e2_08 = fraction( -1674902723.0, 288716400.0); + e2_09 = fraction( -74684743568175.0, 176692375811392.0); + e2_10 = fraction( -734375.0, 4826304.0); + e2_11 = fraction( 171414593.0, 851261400.0); + e2_12 = fraction( 69869.0, 3084480.0); + } + + /** {@inheritDoc} */ + @Override + protected T[] getC() { + + final T sqrt6 = getField().getOne().multiply(6).sqrt(); + + final T[] c = MathArrays.buildArray(getField(), 12); + c[ 0] = sqrt6.add(-6).divide(-67.5); + c[ 1] = sqrt6.add(-6).divide(-45.0); + c[ 2] = sqrt6.add(-6).divide(-30.0); + c[ 3] = sqrt6.add( 6).divide( 30.0); + c[ 4] = fraction(1, 3); + c[ 5] = fraction(1, 4); + c[ 6] = fraction(4, 13); + c[ 7] = fraction(127, 195); + c[ 8] = fraction(3, 5); + c[ 9] = fraction(6, 7); + c[10] = getField().getOne(); + c[11] = getField().getOne(); + + return c; + + } + + /** {@inheritDoc} */ + @Override + protected T[][] getA() { + + final T sqrt6 = getField().getOne().multiply(6).sqrt(); + + final T[][] a = MathArrays.buildArray(getField(), 12, -1); + for (int i = 0; i < a.length; ++i) { + a[i] = MathArrays.buildArray(getField(), i + 1); + } + + a[ 0][ 0] = sqrt6.add(-6).divide(-67.5); + + a[ 1][ 0] = sqrt6.add(-6).divide(-180); + a[ 1][ 1] = sqrt6.add(-6).divide( -40); + + a[ 2][ 0] = sqrt6.add(-6).divide(-120); + a[ 2][ 1] = getField().getZero(); + a[ 2][ 2] = sqrt6.add(-6).divide( -40); + + a[ 3][ 0] = sqrt6.multiply(107).add(462).divide( 3000); + a[ 3][ 1] = getField().getZero(); + a[ 3][ 2] = sqrt6.multiply(197).add(402).divide(-1000); + a[ 3][ 3] = sqrt6.multiply( 73).add(168).divide( 375); + + a[ 4][ 0] = fraction(1, 27); + a[ 4][ 1] = getField().getZero(); + a[ 4][ 2] = getField().getZero(); + a[ 4][ 3] = sqrt6.add( 16).divide( 108); + a[ 4][ 4] = sqrt6.add(-16).divide(-108); + + a[ 5][ 0] = fraction(19, 512); + a[ 5][ 1] = getField().getZero(); + a[ 5][ 2] = getField().getZero(); + a[ 5][ 3] = sqrt6.multiply( 23).add(118).divide(1024); + a[ 5][ 4] = sqrt6.multiply(-23).add(118).divide(1024); + a[ 5][ 5] = fraction(-9, 512); + + a[ 6][ 0] = fraction(13772, 371293); + a[ 6][ 1] = getField().getZero(); + a[ 6][ 2] = getField().getZero(); + a[ 6][ 3] = sqrt6.multiply( 4784).add(51544).divide(371293); + a[ 6][ 4] = sqrt6.multiply(-4784).add(51544).divide(371293); + a[ 6][ 5] = fraction(-5688, 371283); + a[ 6][ 6] = fraction( 3072, 371293); + + a[ 7][ 0] = fraction(58656157643.0, 93983540625.0); + a[ 7][ 1] = getField().getZero(); + a[ 7][ 2] = getField().getZero(); + a[ 7][ 3] = sqrt6.multiply(-318801444819.0).add(-1324889724104.0).divide(626556937500.0); + a[ 7][ 4] = sqrt6.multiply( 318801444819.0).add(-1324889724104.0).divide(626556937500.0); + a[ 7][ 5] = fraction(96044563816.0, 3480871875.0); + a[ 7][ 6] = fraction(5682451879168.0, 281950621875.0); + a[ 7][ 7] = fraction(-165125654.0, 3796875.0); + + a[ 8][ 0] = fraction(8909899.0, 18653125.0); + a[ 8][ 1] = getField().getZero(); + a[ 8][ 2] = getField().getZero(); + a[ 8][ 3] = sqrt6.multiply(-1137963.0).add(-4521408.0).divide(2937500.0); + a[ 8][ 4] = sqrt6.multiply( 1137963.0).add(-4521408.0).divide(2937500.0); + a[ 8][ 5] = fraction(96663078.0, 4553125.0); + a[ 8][ 6] = fraction(2107245056.0, 137915625.0); + a[ 8][ 7] = fraction(-4913652016.0, 147609375.0); + a[ 8][ 8] = fraction(-78894270.0, 3880452869.0); + + a[ 9][ 0] = fraction(-20401265806.0, 21769653311.0); + a[ 9][ 1] = getField().getZero(); + a[ 9][ 2] = getField().getZero(); + a[ 9][ 3] = sqrt6.multiply( 94326.0).add(354216.0).divide(112847.0); + a[ 9][ 4] = sqrt6.multiply(-94326.0).add(354216.0).divide(112847.0); + a[ 9][ 5] = fraction(-43306765128.0, 5313852383.0); + a[ 9][ 6] = fraction(-20866708358144.0, 1126708119789.0); + a[ 9][ 7] = fraction(14886003438020.0, 654632330667.0); + a[ 9][ 8] = fraction(35290686222309375.0, 14152473387134411.0); + a[ 9][ 9] = fraction(-1477884375.0, 485066827.0); + + a[10][ 0] = fraction(39815761.0, 17514443.0); + a[10][ 1] = getField().getZero(); + a[10][ 2] = getField().getZero(); + a[10][ 3] = sqrt6.multiply(-960905.0).add(-3457480.0).divide(551636.0); + a[10][ 4] = sqrt6.multiply( 960905.0).add(-3457480.0).divide(551636.0); + a[10][ 5] = fraction(-844554132.0, 47026969.0); + a[10][ 6] = fraction(8444996352.0, 302158619.0); + a[10][ 7] = fraction(-2509602342.0, 877790785.0); + a[10][ 8] = fraction(-28388795297996250.0, 3199510091356783.0); + 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 + a[11][ 0] = fraction(104257.0, 1920240.0); + a[11][ 1] = getField().getZero(); + a[11][ 2] = getField().getZero(); + a[11][ 3] = getField().getZero(); + a[11][ 4] = getField().getZero(); + a[11][ 5] = fraction(3399327.0, 763840.0); + a[11][ 6] = fraction(66578432.0, 35198415.0); + a[11][ 7] = fraction(-1674902723.0, 288716400.0); + a[11][ 8] = fraction(54980371265625.0, 176692375811392.0); + a[11][ 9] = fraction(-734375.0, 4826304.0); + a[11][10] = fraction(171414593.0, 851261400.0); + a[11][11] = fraction(137909.0, 3084480.0); + + return a; + + } + + /** {@inheritDoc} */ + @Override + protected T[] getB() { + final T[] b = MathArrays.buildArray(getField(), 13); + b[ 0] = fraction(104257, 1929240); + b[ 1] = getField().getZero(); + b[ 2] = getField().getZero(); + b[ 3] = getField().getZero(); + b[ 4] = getField().getZero(); + b[ 5] = fraction( 3399327.0, 763840.0); + b[ 6] = fraction( 66578432.0, 35198415.0); + b[ 7] = fraction( -1674902723.0, 288716400.0); + b[ 8] = fraction( 54980371265625.0, 176692375811392.0); + b[ 9] = fraction( -734375.0, 4826304.0); + b[10] = fraction( 171414593.0, 851261400.0); + b[11] = fraction( 137909.0, 3084480.0); + b[12] = 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); } /** {@inheritDoc} */ @@ -262,22 +364,22 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>> T error2 = h.getField().getZero(); for (int j = 0; j < mainSetDimension; ++j) { - final T errSum1 = yDotK[ 0][j].multiply(E1_01). - add(yDotK[ 5][j].multiply(E1_06)). - add(yDotK[ 6][j].multiply(E1_07)). - add(yDotK[ 7][j].multiply(E1_08)). - add(yDotK[ 8][j].multiply(E1_09)). - add(yDotK[ 9][j].multiply(E1_10)). - add(yDotK[10][j].multiply(E1_11)). - add(yDotK[11][j].multiply(E1_12)); - final T errSum2 = yDotK[ 0][j].multiply(E2_01). - add(yDotK[ 5][j].multiply(E2_06)). - add(yDotK[ 6][j].multiply(E2_07)). - add(yDotK[ 7][j].multiply(E2_08)). - add(yDotK[ 8][j].multiply(E2_09)). - add(yDotK[ 9][j].multiply(E2_10)). - add(yDotK[10][j].multiply(E2_11)). - add(yDotK[11][j].multiply(E2_12)); + final T errSum1 = yDotK[ 0][j].multiply(e1_01). + add(yDotK[ 5][j].multiply(e1_06)). + add(yDotK[ 6][j].multiply(e1_07)). + add(yDotK[ 7][j].multiply(e1_08)). + add(yDotK[ 8][j].multiply(e1_09)). + add(yDotK[ 9][j].multiply(e1_10)). + add(yDotK[10][j].multiply(e1_11)). + add(yDotK[11][j].multiply(e1_12)); + final T errSum2 = yDotK[ 0][j].multiply(e2_01). + add(yDotK[ 5][j].multiply(e2_06)). + add(yDotK[ 6][j].multiply(e2_07)). + add(yDotK[ 7][j].multiply(e2_08)). + add(yDotK[ 8][j].multiply(e2_09)). + add(yDotK[ 9][j].multiply(e2_10)). + add(yDotK[10][j].multiply(e2_11)). + add(yDotK[11][j].multiply(e2_12)); final T yScale = MathUtils.max(y0[j].abs(), y1[j].abs()); final T tol = vecAbsoluteTolerance == null ? http://git-wip-us.apache.org/repos/asf/commons-math/blob/b051dbda/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java index e3a1829..092393a 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java @@ -19,6 +19,9 @@ package org.apache.commons.math3.ode.nonstiff; import org.apache.commons.math3.Field; import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.util.MathArrays; import org.apache.commons.math3.util.MathUtils; @@ -42,30 +45,8 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> /** Integrator method name. */ private static final String METHOD_NAME = "Higham-Hall 5(4)"; - /** Time steps Butcher array. */ - private static final double[] STATIC_C = { - 2.0/9.0, 1.0/3.0, 1.0/2.0, 3.0/5.0, 1.0, 1.0 - }; - - /** Internal weights Butcher array. */ - private static final double[][] STATIC_A = { - {2.0/9.0}, - {1.0/12.0, 1.0/4.0}, - {1.0/8.0, 0.0, 3.0/8.0}, - {91.0/500.0, -27.0/100.0, 78.0/125.0, 8.0/125.0}, - {-11.0/20.0, 27.0/20.0, 12.0/5.0, -36.0/5.0, 5.0}, - {1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0} - }; - - /** Propagation weights Butcher array. */ - private static final double[] STATIC_B = { - 1.0/12.0, 0.0, 27.0/32.0, -4.0/3.0, 125.0/96.0, 5.0/48.0, 0.0 - }; - /** Error weights Butcher array. */ - private static final double[] STATIC_E = { - -1.0/20.0, 0.0, 81.0/160.0, -6.0/5.0, 25.0/32.0, 1.0/16.0, -1.0/10.0 - }; + private final T[] e ; /** Simple constructor. * Build a fifth order Higham and Hall integrator with the given step bounds @@ -83,9 +64,16 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) { - super(field, METHOD_NAME, false, STATIC_C, STATIC_A, STATIC_B, - new HighamHall54FieldStepInterpolator<T>(), + super(field, METHOD_NAME, false, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + e = MathArrays.buildArray(field, 7); + e[0] = fraction(-1, 20); + e[1] = field.getZero(); + e[2] = fraction(81, 160); + e[3] = fraction(-6, 5); + e[4] = fraction(25, 32); + e[5] = fraction( 1, 16); + e[6] = fraction(-1, 10); } /** Simple constructor. @@ -104,9 +92,83 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance) { - super(field, METHOD_NAME, false, STATIC_C, STATIC_A, STATIC_B, - new HighamHall54FieldStepInterpolator<T>(), + super(field, METHOD_NAME, false, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + e = MathArrays.buildArray(field, 7); + e[0] = fraction(-1, 20); + e[1] = field.getZero(); + e[2] = fraction(81, 160); + e[3] = fraction(-6, 5); + e[4] = fraction(25, 32); + e[5] = fraction( 1, 16); + e[6] = fraction(-1, 10); + } + + /** {@inheritDoc} */ + @Override + protected T[] getC() { + final T[] c = MathArrays.buildArray(getField(), 6); + c[0] = fraction(2, 9); + c[1] = fraction(1, 3); + c[2] = fraction(1, 2); + c[3] = fraction(3, 5); + c[4] = getField().getOne(); + c[5] = getField().getOne(); + return c; + } + + /** {@inheritDoc} */ + @Override + protected T[][] getA() { + final T[][] a = MathArrays.buildArray(getField(), 6, -1); + for (int i = 0; i < a.length; ++i) { + a[i] = MathArrays.buildArray(getField(), i + 1); + } + a[0][0] = fraction( 2, 9); + a[1][0] = fraction( 1, 12); + a[1][1] = fraction( 1, 4); + a[2][0] = fraction( 1, 8); + a[2][1] = getField().getZero(); + a[2][2] = fraction( 3, 8); + a[3][0] = fraction( 91, 500); + a[3][1] = fraction( -27, 100); + a[3][2] = fraction( 78, 125); + a[3][3] = fraction( 8, 125); + a[4][0] = fraction( -11, 20); + a[4][1] = fraction( 27, 20); + a[4][2] = fraction( 12, 5); + a[4][3] = fraction( -36, 5); + a[4][4] = fraction( 5, 1); + a[5][0] = fraction( 1, 12); + a[5][1] = getField().getZero(); + a[5][2] = fraction( 27, 32); + a[5][3] = fraction( -4, 3); + a[5][4] = fraction( 125, 96); + a[5][5] = fraction( 5, 48); + return a; + } + + /** {@inheritDoc} */ + @Override + protected T[] getB() { + final T[] b = MathArrays.buildArray(getField(), 7); + b[0] = fraction( 1, 12); + b[1] = getField().getZero(); + b[2] = fraction( 27, 32); + b[3] = fraction( -4, 3); + b[4] = fraction(125, 96); + b[5] = fraction( 5, 48); + b[6] = getField().getZero(); + return b; + } + + /** {@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); } /** {@inheritDoc} */ @@ -122,9 +184,9 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>> T error = getField().getZero(); for (int j = 0; j < mainSetDimension; ++j) { - T errSum = yDotK[0][j].multiply(STATIC_E[0]); - for (int l = 1; l < STATIC_E.length; ++l) { - errSum = errSum.add(yDotK[l][j].multiply(STATIC_E[l])); + T errSum = yDotK[0][j].multiply(e[0]); + for (int l = 1; l < e.length; ++l) { + errSum = errSum.add(yDotK[l][j].multiply(e[l])); } final T yScale = MathUtils.max(y0[j].abs(), y1[j].abs());
