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 5(4) method for solving ODE.
Date Tue, 17 Nov 2015 19:54:32 GMT
Repository: commons-math
Updated Branches:
  refs/heads/field-ode b990f6f2a -> 5bdce0892


Field-based version of Dormand-Prince 5(4) 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/5bdce089
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/5bdce089
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/5bdce089

Branch: refs/heads/field-ode
Commit: 5bdce0892ae4bc516ca92cc2190fb77b84405bdc
Parents: b990f6f
Author: Luc Maisonobe <luc@apache.org>
Authored: Tue Nov 17 20:54:13 2015 +0100
Committer: Luc Maisonobe <luc@apache.org>
Committed: Tue Nov 17 20:54:13 2015 +0100

----------------------------------------------------------------------
 .../DormandPrince54FieldIntegrator.java         | 172 +++++++++++++
 .../DormandPrince54FieldStepInterpolator.java   | 245 +++++++++++++++++++
 2 files changed, 417 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/5bdce089/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
new file mode 100644
index 0000000..7af0375
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
@@ -0,0 +1,172 @@
+/*
+ * 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.MathUtils;
+
+
+/**
+ * This class implements the 5(4) Dormand-Prince integrator for Ordinary
+ * Differential Equations.
+
+ * <p>This integrator is an embedded Runge-Kutta integrator
+ * of order 5(4) 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 7 functions evaluations per step. However, since this
+ * is an <i>fsal</i>, the last evaluation of one step is the same as
+ * the first evaluation of the next step and hence can be avoided. So
+ * the cost is really 6 functions evaluations per step.</p>
+ *
+ * <p>This method has been published (whithout the continuous output
+ * that was added by Shampine in 1986) in the following article :
+ * <pre>
+ *  A family of embedded Runge-Kutta formulae
+ *  J. R. Dormand and P. J. Prince
+ *  Journal of Computational and Applied Mathematics
+ *  volume 6, no 1, 1980, pp. 19-26
+ * </pre></p>
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
+    extends EmbeddedRungeKuttaFieldIntegrator<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;
+
+    // element 2 is zero, so it is neither stored nor used
+
+    /** Error array, element 3. */
+    private static final double E3 =    -71.0 / 16695.0;
+
+    /** Error array, element 4. */
+    private static final double E4 =     71.0 / 1920.0;
+
+    /** Error array, element 5. */
+    private static final double E5 = -17253.0 / 339200.0;
+
+    /** Error array, element 6. */
+    private static final double E6 =     22.0 / 525.0;
+
+    /** Error array, element 7. */
+    private static final double E7 =     -1.0 / 40.0;
+
+    /** Simple constructor.
+     * Build a fifth 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 DormandPrince54FieldIntegrator(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 DormandPrince54FieldStepInterpolator<T>(),
+              minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
+    }
+
+    /** Simple constructor.
+     * Build a fifth 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 DormandPrince54FieldIntegrator(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 DormandPrince54FieldStepInterpolator<T>(),
+              minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int getOrder() {
+        return 5;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
+
+        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 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 ratio  = h.multiply(errSum).divide(tol);
+            error = error.add(ratio.multiply(ratio));
+
+        }
+
+        return error.divide(mainSetDimension).sqrt();
+
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/5bdce089/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
new file mode 100644
index 0000000..1754882
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
@@ -0,0 +1,245 @@
+/*
+ * 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.RealFieldElement;
+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 5(4) Dormand-Prince integrator.
+ *
+ * @see DormandPrince54Integrator
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
+      extends RungeKuttaFieldStepInterpolator<T> {
+
+    /** Last row of the Butcher-array internal weights, element 0. */
+    private static final double A70 =    35.0 /  384.0;
+
+    // element 1 is zero, so it is neither stored nor used
+
+    /** Last row of the Butcher-array internal weights, element 2. */
+    private static final double A72 =   500.0 / 1113.0;
+
+    /** Last row of the Butcher-array internal weights, element 3. */
+    private static final double A73 =   125.0 /  192.0;
+
+    /** Last row of the Butcher-array internal weights, element 4. */
+    private static final double A74 = -2187.0 / 6784.0;
+
+    /** Last row of the Butcher-array internal weights, element 5. */
+    private static final double A75 =    11.0 /   84.0;
+
+    /** Shampine (1986) Dense output, element 0. */
+    private static final double D0 =  -12715105075.0 /  11282082432.0;
+
+    // element 1 is zero, so it is neither stored nor used
+
+    /** Shampine (1986) Dense output, element 2. */
+    private static final double D2 =   87487479700.0 /  32700410799.0;
+
+    /** Shampine (1986) Dense output, element 3. */
+    private static final double D3 =  -10690763975.0 /   1880347072.0;
+
+    /** Shampine (1986) Dense output, element 4. */
+    private static final double D4 =  701980252875.0 / 199316789632.0;
+
+    /** Shampine (1986) Dense output, element 5. */
+    private static final double D5 =   -1453857185.0 /    822651844.0;
+
+    /** Shampine (1986) Dense output, element 6. */
+    private static final double D6 =      69997945.0 /     29380423.0;
+
+    /** First vector for interpolation. */
+    private T[] v1;
+
+    /** Second vector for interpolation. */
+    private T[] v2;
+
+    /** Third vector for interpolation. */
+    private T[] v3;
+
+    /** Fourth vector for interpolation. */
+    private T[] v4;
+
+    /** 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.
+     */
+    DormandPrince54FieldStepInterpolator() {
+        super();
+        v1 = null;
+        v2 = null;
+        v3 = null;
+        v4 = 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
+     */
+    DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T>
interpolator) {
+
+        super(interpolator);
+
+        if (interpolator.v1 == null) {
+
+            v1 = null;
+            v2 = null;
+            v3 = null;
+            v4 = null;
+            vectorsInitialized = false;
+
+        } else {
+
+            v1 = interpolator.v1.clone();
+            v2 = interpolator.v2.clone();
+            v3 = interpolator.v3.clone();
+            v4 = interpolator.v4.clone();
+            vectorsInitialized = interpolator.vectorsInitialized;
+
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected DormandPrince54FieldStepInterpolator<T> doCopy() {
+        return new DormandPrince54FieldStepInterpolator<T>(this);
+    }
+
+
+    /** {@inheritDoc} */
+    @Override
+    protected void reinitialize(final T[] y, final boolean isForward, final FieldEquationsMapper<T>
equationsMapper) {
+        super.reinitialize(y, isForward, equationsMapper);
+        v1 = null;
+        v2 = null;
+        v3 = null;
+        v4 = null;
+        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) {
+
+        if (! vectorsInitialized) {
+
+            if (v1 == null) {
+                v1 = MathArrays.buildArray(time.getField(), previousState.length);
+                v2 = MathArrays.buildArray(time.getField(), previousState.length);
+                v3 = MathArrays.buildArray(time.getField(), previousState.length);
+                v4 = MathArrays.buildArray(time.getField(), previousState.length);
+            }
+
+            // no step finalization is needed for this interpolator
+
+            // we need to compute the interpolation vectors for this time step
+            for (int i = 0; i < previousState.length; ++i) {
+                final T yDot0 = yDotK[0][i];
+                final T yDot2 = yDotK[2][i];
+                final T yDot3 = yDotK[3][i];
+                final T yDot4 = yDotK[4][i];
+                final T yDot5 = yDotK[5][i];
+                final T yDot6 = yDotK[6][i];
+                v1[i] =     yDot0.multiply(A70).
+                        add(yDot2.multiply(A72)).
+                        add(yDot3.multiply(A73)).
+                        add(yDot4.multiply(A74)).
+                        add(yDot5.multiply(A75));
+                v2[i] = yDot0.subtract(v1[i]);
+                v3[i] = v1[i].subtract(v2[i]).subtract(yDot6);
+                v4[i] =     yDot0.multiply(D0).
+                        add(yDot2.multiply(D2)).
+                        add(yDot3.multiply(D3)).
+                        add(yDot4.multiply(D4)).
+                        add(yDot5.multiply(D5)).
+                        add(yDot6.multiply(D6));
+            }
+
+            vectorsInitialized = true;
+
+        }
+
+        // interpolate
+        final T one      = theta.getField().getOne();
+        final T eta      = one.subtract(theta);
+        final T twoTheta = theta.multiply(2);
+        final T dot2     = one.subtract(twoTheta);
+        final T dot3     = theta.multiply(theta.multiply(-3).add(2));
+        final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
+        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(v1[i].
+                                                                     add(eta.multiply(v2[i].
+                                                                                      add(theta.multiply(v3[i].
+                                                                                        
                add(eta.multiply(v4[i])))))))));
+                interpolatedDerivatives[i] =                   v1[i].
+                                             add(dot2.multiply(v2[i])).
+                                             add(dot3.multiply(v3[i])).
+                                             add(dot4.multiply(v4[i]));
+            }
+        } else {
+            for (int i = 0; i < interpolatedState.length; ++i) {
+                interpolatedState[i] = currentState[i].
+                                       subtract(oneMinusThetaH.multiply(v1[i].
+                                                                        subtract(theta.multiply(v2[i].
+                                                                                        
       add(theta.multiply(v3[i].
+                                                                                        
                          add(eta.multiply(v4[i]))))))));
+                interpolatedDerivatives[i] =                   v1[i].
+                                             add(dot2.multiply(v2[i])).
+                                             add(dot3.multiply(v3[i])).
+                                             add(dot4.multiply(v4[i]));
+            }
+        }
+
+        return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+
+    }
+
+}


Mime
View raw message