commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject [math] Field-based version of Dormand-Prince 8(5, 3) method for solving ODE.
Date Sun, 15 Nov 2015 20:58:35 GMT
Repository: commons-math
Updated Branches:
  refs/heads/field-ode e3827069c -> 779e52410


Field-based version of Dormand-Prince 8(5, 3) method for solving ODE.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/779e5241
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/779e5241
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/779e5241

Branch: refs/heads/field-ode
Commit: 779e52410e365f7223bfc5d6e5ad63092a93547f
Parents: e382706
Author: Luc Maisonobe <luc@apache.org>
Authored: Sun Nov 15 21:58:17 2015 +0100
Committer: Luc Maisonobe <luc@apache.org>
Committed: Sun Nov 15 21:58:17 2015 +0100

----------------------------------------------------------------------
 .../DormandPrince853FieldIntegrator.java        | 301 ++++++++++++
 .../DormandPrince853FieldStepInterpolator.java  | 476 +++++++++++++++++++
 2 files changed, 777 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/779e5241/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
new file mode 100644
index 0000000..c4a574c
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
@@ -0,0 +1,301 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+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.util.MathUtils;
+
+
+/**
+ * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary
+ * Differential Equations.
+ *
+ * <p>This integrator is an embedded Runge-Kutta integrator
+ * of order 8(5,3) used in local extrapolation mode (i.e. the solution
+ * is computed using the high order formula) with stepsize control
+ * (and automatic step initialization) and continuous output. This
+ * method uses 12 functions evaluations per step for integration and 4
+ * evaluations for interpolation. However, since the first
+ * interpolation evaluation is the same as the first integration
+ * evaluation of the next step, we have included it in the integrator
+ * rather than in the interpolator and specified the method was an
+ * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
+ * really 12 evaluations per step even if no interpolation is done,
+ * and the overcost of interpolation is only 3 evaluations.</p>
+ *
+ * <p>This method is based on an 8(6) method by Dormand and Prince
+ * (i.e. order 8 for the integration and order 6 for error estimation)
+ * modified by Hairer and Wanner to use a 5th order error estimator
+ * with 3rd order correction. This modification was introduced because
+ * the original method failed in some cases (wrong steps can be
+ * accepted when step size is too large, for example in the
+ * Brusselator problem) and also had <i>severe difficulties when
+ * applied to problems with discontinuities</i>. This modification is
+ * explained in the second edition of the first volume (Nonstiff
+ * Problems) of the reference book by Hairer, Norsett and Wanner:
+ * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
+ * ISBN 3-540-56670-8).</p>
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
+    extends EmbeddedRungeKuttaFieldIntegrator<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;
+
+    // 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;
+
+    /** First error weights array, element 7. */
+    private static final double E1_07 =         -69799717.0 / 140793660.0;
+
+    /** First error weights array, element 8. */
+    private static final double E1_08 =     1230164450203.0 / 739113984000.0;
+
+    /** First error weights array, element 9. */
+    private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;
+
+    /** First error weights array, element 10. */
+    private static final double E1_10 =         464500805.0 / 1389975552.0;
+
+    /** First error weights array, element 11. */
+    private static final double E1_11 =     1606764981773.0 / 19613062656000.0;
+
+    /** First error weights array, element 12. */
+    private static final double E1_12 =           -137909.0 / 6168960.0;
+
+
+    /** Second error weights array, element 1. */
+    private static final double E2_01 =           -364463.0 / 1920240.0;
+
+    // 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;
+
+    /** Second error weights array, element 7. */
+    private static final double E2_07 =          66578432.0 / 35198415.0;
+
+    /** Second error weights array, element 8. */
+    private static final double E2_08 =       -1674902723.0 / 288716400.0;
+
+    /** Second error weights array, element 9. */
+    private static final double E2_09 =   -74684743568175.0 / 176692375811392.0;
+
+    /** Second error weights array, element 10. */
+    private static final double E2_10 =           -734375.0 / 4826304.0;
+
+    /** Second error weights array, element 11. */
+    private static final double E2_11 =         171414593.0 / 851261400.0;
+
+    /** Second error weights array, element 12. */
+    private static final double E2_12 =             69869.0 / 3084480.0;
+
+    /** Simple constructor.
+     * Build an eighth order Dormand-Prince integrator with the given step bounds
+     * @param field field to which the time and state vector elements belong
+     * @param minStep minimal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param maxStep maximal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param scalAbsoluteTolerance allowed absolute error
+     * @param scalRelativeTolerance allowed relative error
+     */
+    public DormandPrince853FieldIntegrator(final Field<T> field,
+                                           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>(),
+              minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+    }
+
+    /** Simple constructor.
+     * Build an eighth order Dormand-Prince integrator with the given step bounds
+     * @param field field to which the time and state vector elements belong
+     * @param minStep minimal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param maxStep maximal step (sign is irrelevant, regardless of
+     * integration direction, forward or backward), the last step can
+     * be smaller than this
+     * @param vecAbsoluteTolerance allowed absolute error
+     * @param vecRelativeTolerance allowed relative error
+     */
+    public DormandPrince853FieldIntegrator(final Field<T> field,
+                                           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>(),
+              minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getOrder() {
+        return 8;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
+        T error1 = h.getField().getZero();
+        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 yScale = MathUtils.max(y0[j].abs(), y1[j].abs());
+            final T tol = vecAbsoluteTolerance == null ?
+                          yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance)
:
+                          yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]);
+            final T ratio1  = errSum1.divide(tol);
+            error1        = error1.add(ratio1.multiply(ratio1));
+            final T ratio2  = errSum2.divide(tol);
+            error2        = error2.add(ratio2.multiply(ratio2));
+        }
+
+        T den = error1.add(error2.multiply(0.01));
+        if (den.getReal() <= 0.0) {
+            den = h.getField().getOne();
+        }
+
+        return h.abs().multiply(error1).divide(den.multiply(mainSetDimension).sqrt());
+
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/779e5241/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
new file mode 100644
index 0000000..b777b4d
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
@@ -0,0 +1,476 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.ode.FieldEquationsMapper;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.util.MathArrays;
+
+/**
+ * This class represents an interpolator over the last step during an
+ * ODE integration for the 8(5,3) Dormand-Prince integrator.
+ *
+ * @see DormandPrince853FieldIntegrator
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
+    extends RungeKuttaFieldStepInterpolator<T> {
+
+    /** Propagation weights, element 1. */
+    private static final double B_01 =         104257.0 / 1920240.0;
+
+    // elements 2 to 5 are zero, so they are neither stored nor used
+
+    /** Propagation weights, element 6. */
+    private static final double B_06 =        3399327.0 / 763840.0;
+
+    /** Propagation weights, element 7. */
+    private static final double B_07 =       66578432.0 / 35198415.0;
+
+    /** Propagation weights, element 8. */
+    private static final double B_08 =    -1674902723.0 / 288716400.0;
+
+    /** Propagation weights, element 9. */
+    private static final double B_09 = 54980371265625.0 / 176692375811392.0;
+
+    /** Propagation weights, element 10. */
+    private static final double B_10 =        -734375.0 / 4826304.0;
+
+    /** Propagation weights, element 11. */
+    private static final double B_11 =      171414593.0 / 851261400.0;
+
+    /** Propagation weights, element 12. */
+    private static final double B_12 =         137909.0 / 3084480.0;
+
+    /** Time step for stage 14 (interpolation only). */
+    private static final double C14    = 1.0 / 10.0;
+
+    /** Internal weights for stage 14, element 1. */
+    private static final double K14_01 =       13481885573.0 / 240030000000.0      - B_01;
+
+    // elements 2 to 5 are zero, so they are neither stored nor used
+
+    /** Internal weights for stage 14, element 6. */
+    private static final double K14_06 =                 0.0                       - B_06;
+
+    /** Internal weights for stage 14, element 7. */
+    private static final double K14_07 =      139418837528.0 / 549975234375.0      - B_07;
+
+    /** Internal weights for stage 14, element 8. */
+    private static final double K14_08 =   -11108320068443.0 / 45111937500000.0    - B_08;
+
+    /** Internal weights for stage 14, element 9. */
+    private static final double K14_09 = -1769651421925959.0 / 14249385146080000.0 - B_09;
+
+    /** Internal weights for stage 14, element 10. */
+    private static final double K14_10 =          57799439.0 / 377055000.0         - B_10;
+
+    /** Internal weights for stage 14, element 11. */
+    private static final double K14_11 =      793322643029.0 / 96734250000000.0    - B_11;
+
+    /** Internal weights for stage 14, element 12. */
+    private static final double K14_12 =        1458939311.0 / 192780000000.0      - B_12;
+
+    /** Internal weights for stage 14, element 13. */
+    private static final double K14_13 =             -4149.0 / 500000.0;
+
+    /** Time step for stage 15 (interpolation only). */
+    private static final double C15    = 1.0 / 5.0;
+
+
+    /** Internal weights for stage 15, element 1. */
+    private static final double K15_01 =     1595561272731.0 / 50120273500000.0    - B_01;
+
+    // elements 2 to 5 are zero, so they are neither stored nor used
+
+    /** Internal weights for stage 15, element 6. */
+    private static final double K15_06 =      975183916491.0 / 34457688031250.0    - B_06;
+
+    /** Internal weights for stage 15, element 7. */
+    private static final double K15_07 =    38492013932672.0 / 718912673015625.0   - B_07;
+
+    /** Internal weights for stage 15, element 8. */
+    private static final double K15_08 = -1114881286517557.0 / 20298710767500000.0 - B_08;
+
+    /** Internal weights for stage 15, element 9. */
+    private static final double K15_09 =                 0.0                       - B_09;
+
+    /** Internal weights for stage 15, element 10. */
+    private static final double K15_10 =                 0.0                       - B_10;
+
+    /** Internal weights for stage 15, element 11. */
+    private static final double K15_11 =    -2538710946863.0 / 23431227861250000.0 - B_11;
+
+    /** Internal weights for stage 15, element 12. */
+    private static final double K15_12 =        8824659001.0 / 23066716781250.0    - B_12;
+
+    /** Internal weights for stage 15, element 13. */
+    private static final double K15_13 =      -11518334563.0 / 33831184612500.0;
+
+    /** Internal weights for stage 15, element 14. */
+    private static final double K15_14 =        1912306948.0 / 13532473845.0;
+
+    /** Time step for stage 16 (interpolation only). */
+    private static final double C16    = 7.0 / 9.0;
+
+
+    /** Internal weights for stage 16, element 1. */
+    private static final double K16_01 =      -13613986967.0 / 31741908048.0       - B_01;
+
+    // elements 2 to 5 are zero, so they are neither stored nor used
+
+    /** Internal weights for stage 16, element 6. */
+    private static final double K16_06 =       -4755612631.0 / 1012344804.0        - B_06;
+
+    /** Internal weights for stage 16, element 7. */
+    private static final double K16_07 =    42939257944576.0 / 5588559685701.0     - B_07;
+
+    /** Internal weights for stage 16, element 8. */
+    private static final double K16_08 =    77881972900277.0 / 19140370552944.0    - B_08;
+
+    /** Internal weights for stage 16, element 9. */
+    private static final double K16_09 =    22719829234375.0 / 63689648654052.0    - B_09;
+
+    /** Internal weights for stage 16, element 10. */
+    private static final double K16_10 =                 0.0                       - B_10;
+
+    /** Internal weights for stage 16, element 11. */
+    private static final double K16_11 =                 0.0                       - B_11;
+
+    /** Internal weights for stage 16, element 12. */
+    private static final double K16_12 =                 0.0                       - B_12;
+
+    /** Internal weights for stage 16, element 13. */
+    private static final double K16_13 =       -1199007803.0 / 857031517296.0;
+
+    /** Internal weights for stage 16, element 14. */
+    private static final double K16_14 =      157882067000.0 / 53564469831.0;
+
+    /** Internal weights for stage 16, element 15. */
+    private static final double K16_15 =     -290468882375.0 / 31741908048.0;
+
+    /** Interpolation weights.
+     * (beware that only the non-null values are in the table)
+     */
+    private static final double[][] D = {
+
+      {        -17751989329.0 / 2106076560.0,               4272954039.0 / 7539864640.0,
+              -118476319744.0 / 38604839385.0,            755123450731.0 / 316657731600.0,
+        3692384461234828125.0 / 1744130441634250432.0,     -4612609375.0 / 5293382976.0,
+              2091772278379.0 / 933644586600.0,             2136624137.0 / 3382989120.0,
+                    -126493.0 / 1421424.0,                    98350000.0 / 5419179.0,
+                  -18878125.0 / 2053168.0,                 -1944542619.0 / 438351368.0},
+
+      {         32941697297.0 / 3159114840.0,             456696183123.0 / 1884966160.0,
+             19132610714624.0 / 115814518155.0,       -177904688592943.0 / 474986597400.0,
+       -4821139941836765625.0 / 218016305204281304.0,      30702015625.0 / 3970037232.0,
+            -85916079474274.0 / 2800933759800.0,           -5919468007.0 / 634310460.0,
+                    2479159.0 / 157936.0,                    -18750000.0 / 602131.0,
+                  -19203125.0 / 2053168.0,                 15700361463.0 / 438351368.0},
+
+      {         12627015655.0 / 631822968.0,              -72955222965.0 / 188496616.0,
+            -13145744952320.0 / 69488710893.0,          30084216194513.0 / 56998391688.0,
+        -296858761006640625.0 / 25648977082856624.0,         569140625.0 / 82709109.0,
+               -18684190637.0 / 18672891732.0,                69644045.0 / 89549712.0,
+                  -11847025.0 / 4264272.0,                  -978650000.0 / 16257537.0,
+                  519371875.0 / 6159504.0,                  5256837225.0 / 438351368.0},
+
+      {          -450944925.0 / 17550638.0,               -14532122925.0 / 94248308.0,
+              -595876966400.0 / 2573655959.0,             188748653015.0 / 527762886.0,
+        2545485458115234375.0 / 27252038150535163.0,       -1376953125.0 / 36759604.0,
+                53995596795.0 / 518691437.0,                 210311225.0 / 7047894.0,
+                   -1718875.0 / 39484.0,                      58000000.0 / 602131.0,
+                   -1546875.0 / 39484.0,                   -1262172375.0 / 8429834.0}
+
+    };
+
+    /** Last evaluations. */
+    private T[][] yDotKLast;
+
+    /** Vectors for interpolation. */
+    private T[][] v;
+
+    /** Initialization indicator for the interpolation vectors. */
+    private boolean vectorsInitialized;
+
+    /** Simple constructor.
+     * This constructor builds an instance that is not usable yet, the
+     * {@link #reinitialize} method should be called before using the
+     * instance in order to initialize the internal arrays. This
+     * constructor is used only in order to delay the initialization in
+     * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
+     * prototyping design pattern to create the step interpolators by
+     * cloning an uninitialized model and latter initializing the copy.
+     */
+    DormandPrince853FieldStepInterpolator() {
+        super();
+        yDotKLast = null;
+        v         = null;
+        vectorsInitialized = false;
+    }
+
+    /** Copy constructor.
+     * @param interpolator interpolator to copy from. The copy is a deep
+     * copy: its arrays are separated from the original arrays of the
+     * instance
+     */
+    DormandPrince853FieldStepInterpolator(final DormandPrince853FieldStepInterpolator<T>
interpolator) {
+
+        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;
+
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected DormandPrince853FieldStepInterpolator<T> doCopy() {
+        return new DormandPrince853FieldStepInterpolator<T>(this);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void reinitialize(final T[] y, final boolean isForward, final FieldEquationsMapper<T>
equationsMapper) {
+
+        super.reinitialize(y, isForward, equationsMapper);
+
+        final int dimension = currentState.length;
+
+        yDotKLast = MathArrays.buildArray(y[0].getField(), 3, dimension);
+        v         = MathArrays.buildArray(y[0].getField(), 7, dimension);
+
+        vectorsInitialized = false;
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void storeState(final FieldODEStateAndDerivative<T> state) {
+        super.storeState(state);
+        vectorsInitialized = false;
+    }
+
+    /** {@inheritDoc} */
+    @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);
+        final T theta2   = theta.multiply(theta);
+        final T dot1     = one.subtract(twoTheta);
+        final T dot2     = theta.multiply(theta.multiply(-3).add(2));
+        final T dot3     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
+        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);
+
+        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]));
+            }
+        } 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]));
+            }
+        }
+
+        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);
+
+    }
+
+}


Mime
View raw message