DormandPrince54Integrator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.ode.nonstiff;

  20. import org.apache.commons.math4.core.jdkmath.JdkMath;


  23. /**
  24.  * This class implements the 5(4) Dormand-Prince integrator for Ordinary
  25.  * Differential Equations.

  27.  * <p>This integrator is an embedded Runge-Kutta integrator
  28.  * of order 5(4) used in local extrapolation mode (i.e. the solution
  29.  * is computed using the high order formula) with stepsize control
  30.  * (and automatic step initialization) and continuous output. This
  31.  * method uses 7 functions evaluations per step. However, since this
  32.  * is an <i>fsal</i>, the last evaluation of one step is the same as
  33.  * the first evaluation of the next step and hence can be avoided. So
  34.  * the cost is really 6 functions evaluations per step.</p>
  35.  *
  36.  * <p>This method has been published (whithout the continuous output
  37.  * that was added by Shampine in 1986) in the following article :
  38.  * <pre>
  39.  *  A family of embedded Runge-Kutta formulae
  40.  *  J. R. Dormand and P. J. Prince
  41.  *  Journal of Computational and Applied Mathematics
  42.  *  volume 6, no 1, 1980, pp. 19-26
  43.  * </pre>
  44.  *
  45.  * @since 1.2
  46.  */

  48. public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator {

  50.   /** Integrator method name. */
  51.   private static final String METHOD_NAME = "Dormand-Prince 5(4)";

  53.   /** Time steps Butcher array. */
  54.   private static final double[] STATIC_C = {
  55.     1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0
  56.   };

  58.   /** Internal weights Butcher array. */
  59.   private static final double[][] STATIC_A = {
  60.     {1.0/5.0},
  61.     {3.0/40.0, 9.0/40.0},
  62.     {44.0/45.0, -56.0/15.0, 32.0/9.0},
  63.     {19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0,  -212.0/729.0},
  64.     {9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0},
  65.     {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0}
  66.   };

  68.   /** Propagation weights Butcher array. */
  69.   private static final double[] STATIC_B = {
  70.     35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0
  71.   };

  73.   /** Error array, element 1. */
  74.   private static final double E1 =     71.0 / 57600.0;

  76.   // element 2 is zero, so it is neither stored nor used

  78.   /** Error array, element 3. */
  79.   private static final double E3 =    -71.0 / 16695.0;

  81.   /** Error array, element 4. */
  82.   private static final double E4 =     71.0 / 1920.0;

  84.   /** Error array, element 5. */
  85.   private static final double E5 = -17253.0 / 339200.0;

  87.   /** Error array, element 6. */
  88.   private static final double E6 =     22.0 / 525.0;

  90.   /** Error array, element 7. */
  91.   private static final double E7 =     -1.0 / 40.0;

  93.   /** Simple constructor.
  94.    * Build a fifth order Dormand-Prince integrator with the given step bounds
  95.    * @param minStep minimal step (sign is irrelevant, regardless of
  96.    * integration direction, forward or backward), the last step can
  97.    * be smaller than this
  98.    * @param maxStep maximal step (sign is irrelevant, regardless of
  99.    * integration direction, forward or backward), the last step can
  100.    * be smaller than this
  101.    * @param scalAbsoluteTolerance allowed absolute error
  102.    * @param scalRelativeTolerance allowed relative error
  103.    */
  104.   public DormandPrince54Integrator(final double minStep, final double maxStep,
  105.                                    final double scalAbsoluteTolerance,
  106.                                    final double scalRelativeTolerance) {
  107.     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, new DormandPrince54StepInterpolator(),
  108.           minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  109.   }

  111.   /** Simple constructor.
  112.    * Build a fifth order Dormand-Prince integrator with the given step bounds
  113.    * @param minStep minimal step (sign is irrelevant, regardless of
  114.    * integration direction, forward or backward), the last step can
  115.    * be smaller than this
  116.    * @param maxStep maximal step (sign is irrelevant, regardless of
  117.    * integration direction, forward or backward), the last step can
  118.    * be smaller than this
  119.    * @param vecAbsoluteTolerance allowed absolute error
  120.    * @param vecRelativeTolerance allowed relative error
  121.    */
  122.   public DormandPrince54Integrator(final double minStep, final double maxStep,
  123.                                    final double[] vecAbsoluteTolerance,
  124.                                    final double[] vecRelativeTolerance) {
  125.     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B, new DormandPrince54StepInterpolator(),
  126.           minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  127.   }

  129.   /** {@inheritDoc} */
  130.   @Override
  131.   public int getOrder() {
  132.     return 5;
  133.   }

  135.   /** {@inheritDoc} */
  136.   @Override
  137.   protected double estimateError(final double[][] yDotK,
  138.                                  final double[] y0, final double[] y1,
  139.                                  final double h) {

  141.     double error = 0;

  143.     for (int j = 0; j < mainSetDimension; ++j) {
  144.         final double errSum = E1 * yDotK[0][j] +  E3 * yDotK[2][j] +
  145.                               E4 * yDotK[3][j] +  E5 * yDotK[4][j] +
  146.                               E6 * yDotK[5][j] +  E7 * yDotK[6][j];

  148.         final double yScale = JdkMath.max(JdkMath.abs(y0[j]), JdkMath.abs(y1[j]));
  149.         final double tol = (vecAbsoluteTolerance == null) ?
  150.                            (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
  151.                                (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
  152.         final double ratio  = h * errSum / tol;
  153.         error += ratio * ratio;
  154.     }

  156.     return JdkMath.sqrt(error / mainSetDimension);
  157.   }
  158. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT