DormandPrince853Integrator.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 8(5,3) Dormand-Prince integrator for Ordinary
  25.  * Differential Equations.
  26.  *
  27.  * <p>This integrator is an embedded Runge-Kutta integrator
  28.  * of order 8(5,3) 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 12 functions evaluations per step for integration and 4
  32.  * evaluations for interpolation. However, since the first
  33.  * interpolation evaluation is the same as the first integration
  34.  * evaluation of the next step, we have included it in the integrator
  35.  * rather than in the interpolator and specified the method was an
  36.  * <i>fsal</i>. Hence, despite we have 13 stages here, the cost is
  37.  * really 12 evaluations per step even if no interpolation is done,
  38.  * and the overcost of interpolation is only 3 evaluations.</p>
  39.  *
  40.  * <p>This method is based on an 8(6) method by Dormand and Prince
  41.  * (i.e. order 8 for the integration and order 6 for error estimation)
  42.  * modified by Hairer and Wanner to use a 5th order error estimator
  43.  * with 3rd order correction. This modification was introduced because
  44.  * the original method failed in some cases (wrong steps can be
  45.  * accepted when step size is too large, for example in the
  46.  * Brusselator problem) and also had <i>severe difficulties when
  47.  * applied to problems with discontinuities</i>. This modification is
  48.  * explained in the second edition of the first volume (Nonstiff
  49.  * Problems) of the reference book by Hairer, Norsett and Wanner:
  50.  * <i>Solving Ordinary Differential Equations</i> (Springer-Verlag,
  51.  * ISBN 3-540-56670-8).</p>
  52.  *
  53.  * @since 1.2
  54.  */

  56. public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {

  58.   /** Integrator method name. */
  59.   private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";

  61.   /** Time steps Butcher array. */
  62.   private static final double[] STATIC_C = {
  63.     (12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0, (6.0 - JdkMath.sqrt(6.0)) / 45.0, (6.0 - JdkMath.sqrt(6.0)) / 30.0,
  64.     (6.0 + JdkMath.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,
  65.     6.0/7.0, 1.0, 1.0
  66.   };

  68.   /** Internal weights Butcher array. */
  69.   private static final double[][] STATIC_A = {

  71.     // k2
  72.     {(12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0},

  74.     // k3
  75.     {(6.0 - JdkMath.sqrt(6.0)) / 180.0, (6.0 - JdkMath.sqrt(6.0)) / 60.0},

  77.     // k4
  78.     {(6.0 - JdkMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - JdkMath.sqrt(6.0)) / 40.0},

  80.     // k5
  81.     {(462.0 + 107.0 * JdkMath.sqrt(6.0)) / 3000.0, 0.0,
  82.      (-402.0 - 197.0 * JdkMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * JdkMath.sqrt(6.0)) / 375.0},

  84.     // k6
  85.     {1.0 / 27.0, 0.0, 0.0, (16.0 + JdkMath.sqrt(6.0)) / 108.0, (16.0 - JdkMath.sqrt(6.0)) / 108.0},

  87.     // k7
  88.     {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * JdkMath.sqrt(6.0)) / 1024.0,
  89.      (118.0 - 23.0 * JdkMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0},

  91.     // k8
  92.     {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * JdkMath.sqrt(6.0)) / 371293.0,
  93.      (51544.0 - 4784.0 * JdkMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},

  95.     // k9
  96.     {58656157643.0 / 93983540625.0, 0.0, 0.0,
  97.      (-1324889724104.0 - 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0,
  98.      (-1324889724104.0 + 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0,
  99.      96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,
  100.      -165125654.0 / 3796875.0},

  102.     // k10
  103.     {8909899.0 / 18653125.0, 0.0, 0.0,
  104.      (-4521408.0 - 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0,
  105.      (-4521408.0 + 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0,
  106.      96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,
  107.      -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},

  109.     // k11
  110.     {-20401265806.0 / 21769653311.0, 0.0, 0.0,
  111.      (354216.0 + 94326.0 * JdkMath.sqrt(6.0)) / 112847.0,
  112.      (354216.0 - 94326.0 * JdkMath.sqrt(6.0)) / 112847.0,
  113.      -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,
  114.      14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,
  115.      -1477884375.0 / 485066827.0},

  117.     // k12
  118.     {39815761.0 / 17514443.0, 0.0, 0.0,
  119.      (-3457480.0 - 960905.0 * JdkMath.sqrt(6.0)) / 551636.0,
  120.      (-3457480.0 + 960905.0 * JdkMath.sqrt(6.0)) / 551636.0,
  121.      -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,
  122.      -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,
  123.      226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},

  125.     // k13 should be for interpolation only, but since it is the same
  126.     // stage as the first evaluation of the next step, we perform it
  127.     // here at no cost by specifying this is an fsal method
  128.     {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,
  129.      66578432.0/35198415.0, -1674902723.0/288716400.0,
  130.      54980371265625.0/176692375811392.0, -734375.0/4826304.0,
  131.      171414593.0/851261400.0, 137909.0/3084480.0}
  132.   };

  134.   /** Propagation weights Butcher array. */
  135.   private static final double[] STATIC_B = {
  136.       104257.0/1920240.0,
  137.       0.0,
  138.       0.0,
  139.       0.0,
  140.       0.0,
  141.       3399327.0/763840.0,
  142.       66578432.0/35198415.0,
  143.       -1674902723.0/288716400.0,
  144.       54980371265625.0/176692375811392.0,
  145.       -734375.0/4826304.0,
  146.       171414593.0/851261400.0,
  147.       137909.0/3084480.0,
  148.       0.0
  149.   };

  151.   /** First error weights array, element 1. */
  152.   private static final double E1_01 =         116092271.0 / 8848465920.0;

  154.   // elements 2 to 5 are zero, so they are neither stored nor used

  156.   /** First error weights array, element 6. */
  157.   private static final double E1_06 =          -1871647.0 / 1527680.0;

  159.   /** First error weights array, element 7. */
  160.   private static final double E1_07 =         -69799717.0 / 140793660.0;

  162.   /** First error weights array, element 8. */
  163.   private static final double E1_08 =     1230164450203.0 / 739113984000.0;

  165.   /** First error weights array, element 9. */
  166.   private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;

  168.   /** First error weights array, element 10. */
  169.   private static final double E1_10 =         464500805.0 / 1389975552.0;

  171.   /** First error weights array, element 11. */
  172.   private static final double E1_11 =     1606764981773.0 / 19613062656000.0;

  174.   /** First error weights array, element 12. */
  175.   private static final double E1_12 =           -137909.0 / 6168960.0;


  178.   /** Second error weights array, element 1. */
  179.   private static final double E2_01 =           -364463.0 / 1920240.0;

  181.   // elements 2 to 5 are zero, so they are neither stored nor used

  183.   /** Second error weights array, element 6. */
  184.   private static final double E2_06 =           3399327.0 / 763840.0;

  186.   /** Second error weights array, element 7. */
  187.   private static final double E2_07 =          66578432.0 / 35198415.0;

  189.   /** Second error weights array, element 8. */
  190.   private static final double E2_08 =       -1674902723.0 / 288716400.0;

  192.   /** Second error weights array, element 9. */
  193.   private static final double E2_09 =   -74684743568175.0 / 176692375811392.0;

  195.   /** Second error weights array, element 10. */
  196.   private static final double E2_10 =           -734375.0 / 4826304.0;

  198.   /** Second error weights array, element 11. */
  199.   private static final double E2_11 =         171414593.0 / 851261400.0;

  201.   /** Second error weights array, element 12. */
  202.   private static final double E2_12 =             69869.0 / 3084480.0;

  204.   /** Simple constructor.
  205.    * Build an eighth order Dormand-Prince integrator with the given step bounds
  206.    * @param minStep minimal step (sign is irrelevant, regardless of
  207.    * integration direction, forward or backward), the last step can
  208.    * be smaller than this
  209.    * @param maxStep maximal step (sign is irrelevant, regardless of
  210.    * integration direction, forward or backward), the last step can
  211.    * be smaller than this
  212.    * @param scalAbsoluteTolerance allowed absolute error
  213.    * @param scalRelativeTolerance allowed relative error
  214.    */
  215.   public DormandPrince853Integrator(final double minStep, final double maxStep,
  216.                                     final double scalAbsoluteTolerance,
  217.                                     final double scalRelativeTolerance) {
  218.     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
  219.           new DormandPrince853StepInterpolator(),
  220.           minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  221.   }

  223.   /** Simple constructor.
  224.    * Build an eighth order Dormand-Prince integrator with the given step bounds
  225.    * @param minStep minimal step (sign is irrelevant, regardless of
  226.    * integration direction, forward or backward), the last step can
  227.    * be smaller than this
  228.    * @param maxStep maximal step (sign is irrelevant, regardless of
  229.    * integration direction, forward or backward), the last step can
  230.    * be smaller than this
  231.    * @param vecAbsoluteTolerance allowed absolute error
  232.    * @param vecRelativeTolerance allowed relative error
  233.    */
  234.   public DormandPrince853Integrator(final double minStep, final double maxStep,
  235.                                     final double[] vecAbsoluteTolerance,
  236.                                     final double[] vecRelativeTolerance) {
  237.     super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
  238.           new DormandPrince853StepInterpolator(),
  239.           minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  240.   }

  242.   /** {@inheritDoc} */
  243.   @Override
  244.   public int getOrder() {
  245.     return 8;
  246.   }

  248.   /** {@inheritDoc} */
  249.   @Override
  250.   protected double estimateError(final double[][] yDotK,
  251.                                  final double[] y0, final double[] y1,
  252.                                  final double h) {
  253.     double error1 = 0;
  254.     double error2 = 0;

  256.     for (int j = 0; j < mainSetDimension; ++j) {
  257.       final double errSum1 = E1_01 * yDotK[0][j]  + E1_06 * yDotK[5][j] +
  258.                              E1_07 * yDotK[6][j]  + E1_08 * yDotK[7][j] +
  259.                              E1_09 * yDotK[8][j]  + E1_10 * yDotK[9][j] +
  260.                              E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
  261.       final double errSum2 = E2_01 * yDotK[0][j]  + E2_06 * yDotK[5][j] +
  262.                              E2_07 * yDotK[6][j]  + E2_08 * yDotK[7][j] +
  263.                              E2_09 * yDotK[8][j]  + E2_10 * yDotK[9][j] +
  264.                              E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];

  266.       final double yScale = JdkMath.max(JdkMath.abs(y0[j]), JdkMath.abs(y1[j]));
  267.       final double tol = (vecAbsoluteTolerance == null) ?
  268.                          (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
  269.                          (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
  270.       final double ratio1  = errSum1 / tol;
  271.       error1        += ratio1 * ratio1;
  272.       final double ratio2  = errSum2 / tol;
  273.       error2        += ratio2 * ratio2;
  274.     }

  276.     double den = error1 + 0.01 * error2;
  277.     if (den <= 0.0) {
  278.       den = 1.0;
  279.     }

  281.     return JdkMath.abs(h) * error1 / JdkMath.sqrt(mainSetDimension * den);
  282.   }
  283. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT