DormandPrince853Integrator.java

  1. /*
  2.  * Licensed to the Hipparchus project 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 Hipparchus project 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.  *      https://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.hipparchus.ode.nonstiff;

  20. import org.hipparchus.ode.EquationsMapper;
  21. import org.hipparchus.ode.ODEStateAndDerivative;
  22. import org.hipparchus.ode.nonstiff.interpolators.DormandPrince853StateInterpolator;
  23. import org.hipparchus.util.FastMath;


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

  58. public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {

  60.     /** Name of integration scheme. */
  61.     public static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";

  63.     /** First error weights array, element 1. */
  64.     static final double E1_01 =         116092271.0 / 8848465920.0;

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

  68.     /** First error weights array, element 6. */
  69.     static final double E1_06 =          -1871647.0 / 1527680.0;

  71.     /** First error weights array, element 7. */
  72.     static final double E1_07 =         -69799717.0 / 140793660.0;

  74.     /** First error weights array, element 8. */
  75.     static final double E1_08 =     1230164450203.0 / 739113984000.0;

  77.     /** First error weights array, element 9. */
  78.     static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;

  80.     /** First error weights array, element 10. */
  81.     static final double E1_10 =         464500805.0 / 1389975552.0;

  83.     /** First error weights array, element 11. */
  84.     static final double E1_11 =     1606764981773.0 / 19613062656000.0;

  86.     /** First error weights array, element 12. */
  87.     static final double E1_12 =           -137909.0 / 6168960.0;


  90.     /** Second error weights array, element 1. */
  91.     static final double E2_01 =           -364463.0 / 1920240.0;

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

  95.     /** Second error weights array, element 6. */
  96.     static final double E2_06 =           3399327.0 / 763840.0;

  98.     /** Second error weights array, element 7. */
  99.     static final double E2_07 =          66578432.0 / 35198415.0;

  101.     /** Second error weights array, element 8. */
  102.     static final double E2_08 =       -1674902723.0 / 288716400.0;

  104.     /** Second error weights array, element 9. */
  105.     static final double E2_09 =   -74684743568175.0 / 176692375811392.0;

  107.     /** Second error weights array, element 10. */
  108.     static final double E2_10 =           -734375.0 / 4826304.0;

  110.     /** Second error weights array, element 11. */
  111.     static final double E2_11 =         171414593.0 / 851261400.0;

  113.     /** Second error weights array, element 12. */
  114.     static final double E2_12 =             69869.0 / 3084480.0;

  116.     /** Simple constructor.
  117.      * Build a fifth order Dormand-Prince integrator with the given step bounds
  118.      * @param minStep minimal step (sign is irrelevant, regardless of
  119.      * integration direction, forward or backward), the last step can
  120.      * be smaller than this
  121.      * @param maxStep maximal step (sign is irrelevant, regardless of
  122.      * integration direction, forward or backward), the last step can
  123.      * be smaller than this
  124.      * @param scalAbsoluteTolerance allowed absolute error
  125.      * @param scalRelativeTolerance allowed relative error
  126.      */
  127.     public DormandPrince853Integrator(final double minStep, final double maxStep,
  128.                                       final double scalAbsoluteTolerance,
  129.                                       final double scalRelativeTolerance) {
  130.         super(METHOD_NAME, 12,
  131.               minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  132.     }

  134.     /** Simple constructor.
  135.      * Build a fifth order Dormand-Prince integrator with the given step bounds
  136.      * @param minStep minimal step (sign is irrelevant, regardless of
  137.      * integration direction, forward or backward), the last step can
  138.      * be smaller than this
  139.      * @param maxStep maximal step (sign is irrelevant, regardless of
  140.      * integration direction, forward or backward), the last step can
  141.      * be smaller than this
  142.      * @param vecAbsoluteTolerance allowed absolute error
  143.      * @param vecRelativeTolerance allowed relative error
  144.      */
  145.     public DormandPrince853Integrator(final double minStep, final double maxStep,
  146.                                       final double[] vecAbsoluteTolerance,
  147.                                       final double[] vecRelativeTolerance) {
  148.         super(METHOD_NAME, 12,
  149.               minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  150.     }

  152.     /** {@inheritDoc} */
  153.     @Override
  154.     public double[] getC() {
  155.         final double sqrt6 = FastMath.sqrt(6.0);
  156.         return new double[] {
  157.             (12.0 - 2.0 * sqrt6) / 135.0,
  158.             (6.0 - sqrt6) / 45.0,
  159.             (6.0 - sqrt6) / 30.0,
  160.             (6.0 + sqrt6) / 30.0,
  161.             1.0/3.0,
  162.             1.0/4.0,
  163.             4.0/13.0,
  164.             127.0/195.0,
  165.             3.0/5.0,
  166.             6.0/7.0,
  167.             1.0,
  168.             1.0,
  169.             1.0/10.0,
  170.             1.0/5.0,
  171.             7.0/9.0
  172.         };
  173.     }

  175.     /** {@inheritDoc} */
  176.     @Override
  177.     public double[][] getA() {
  178.         final double sqrt6 = FastMath.sqrt(6.0);
  179.         return new double[][] {
  180.             {
  181.                 (12.0 - 2.0 * sqrt6) / 135.0
  182.             }, {
  183.                 (6.0 - sqrt6) / 180.0,
  184.                 (6.0 - sqrt6) / 60.0
  185.             }, {
  186.                 (6.0 - sqrt6) / 120.0,
  187.                 0.0,
  188.                 (6.0 - sqrt6) / 40.0
  189.             }, {
  190.                 (462.0 + 107.0 * sqrt6) / 3000.0,
  191.                 0.0,
  192.                 (-402.0 - 197.0 * sqrt6) / 1000.0,
  193.                 (168.0 + 73.0 * sqrt6) / 375.0
  194.             }, {
  195.                 1.0 / 27.0,
  196.                 0.0,
  197.                 0.0,
  198.                 (16.0 + sqrt6) / 108.0,
  199.                 (16.0 - sqrt6) / 108.0
  200.             }, {
  201.                 19.0 / 512.0,
  202.                 0.0,
  203.                 0.0,
  204.                 (118.0 + 23.0 * sqrt6) / 1024.0,
  205.                 (118.0 - 23.0 * sqrt6) / 1024.0,
  206.                 -9.0 / 512.0
  207.             }, {
  208.                 13772.0 / 371293.0,
  209.                 0.0,
  210.                 0.0,
  211.                 (51544.0 + 4784.0 * sqrt6) / 371293.0,
  212.                 (51544.0 - 4784.0 * sqrt6) / 371293.0,
  213.                 -5688.0 / 371293.0,
  214.                 3072.0 / 371293.0
  215.             }, {
  216.                 58656157643.0 / 93983540625.0,
  217.                 0.0,
  218.                 0.0,
  219.                 (-1324889724104.0 - 318801444819.0 * sqrt6) / 626556937500.0,
  220.                 (-1324889724104.0 + 318801444819.0 * sqrt6) / 626556937500.0,
  221.                 96044563816.0 / 3480871875.0,
  222.                 5682451879168.0 / 281950621875.0,
  223.                 -165125654.0 / 3796875.0
  224.             }, {
  225.                 8909899.0 / 18653125.0,
  226.                 0.0,
  227.                 0.0,
  228.                 (-4521408.0 - 1137963.0 * sqrt6) / 2937500.0,
  229.                 (-4521408.0 + 1137963.0 * sqrt6) / 2937500.0,
  230.                 96663078.0 / 4553125.0,
  231.                 2107245056.0 / 137915625.0,
  232.                 -4913652016.0 / 147609375.0,
  233.                 -78894270.0 / 3880452869.0
  234.             }, {
  235.                 -20401265806.0 / 21769653311.0,
  236.                 0.0,
  237.                 0.0,
  238.                 (354216.0 + 94326.0 * sqrt6) / 112847.0,
  239.                 (354216.0 - 94326.0 * sqrt6) / 112847.0,
  240.                 -43306765128.0 / 5313852383.0,
  241.                 -20866708358144.0 / 1126708119789.0,
  242.                 14886003438020.0 / 654632330667.0,
  243.                 35290686222309375.0 / 14152473387134411.0,
  244.                 -1477884375.0 / 485066827.0
  245.             }, {
  246.                 39815761.0 / 17514443.0,
  247.                 0.0,
  248.                 0.0,
  249.                 (-3457480.0 - 960905.0 * sqrt6) / 551636.0,
  250.                 (-3457480.0 + 960905.0 * sqrt6) / 551636.0,
  251.                 -844554132.0 / 47026969.0,
  252.                 8444996352.0 / 302158619.0,
  253.                 -2509602342.0 / 877790785.0,
  254.                 -28388795297996250.0 / 3199510091356783.0,
  255.                 226716250.0 / 18341897.0,
  256.                 1371316744.0 / 2131383595.0
  257.             }, {
  258.                 // the following stage is both for interpolation and the first stage in next step
  259.                 // (the coefficients are identical to the B array)
  260.                 104257.0/1920240.0,
  261.                 0.0,
  262.                 0.0,
  263.                 0.0,
  264.                 0.0,
  265.                 3399327.0/763840.0,
  266.                 66578432.0/35198415.0,
  267.                 -1674902723.0/288716400.0,
  268.                 54980371265625.0/176692375811392.0,
  269.                 -734375.0/4826304.0,
  270.                 171414593.0/851261400.0,
  271.                 137909.0/3084480.0
  272.             }, {
  273.                 // the following stages are for interpolation only
  274.                 13481885573.0 / 240030000000.0,
  275.                 0.0,
  276.                 0.0,
  277.                 0.0,
  278.                 0.0,
  279.                 0.0,
  280.                 139418837528.0 / 549975234375.0,
  281.                 -11108320068443.0 / 45111937500000.0,
  282.                 -1769651421925959.0 / 14249385146080000.0,
  283.                 57799439.0 / 377055000.0,
  284.                 793322643029.0 / 96734250000000.0,
  285.                 1458939311.0 / 192780000000.0,
  286.                 -4149.0 / 500000.0
  287.             }, {
  288.                 1595561272731.0 / 50120273500000.0,
  289.                 0.0,
  290.                 0.0,
  291.                 0.0,
  292.                 0.0,
  293.                 975183916491.0 /  34457688031250.0,
  294.                 38492013932672.0 /  718912673015625.0,
  295.                 -1114881286517557.0 /  20298710767500000.0,
  296.                 0.0,
  297.                 0.0,
  298.                 -2538710946863.0 /  23431227861250000.0,
  299.                 8824659001.0 /  23066716781250.0,
  300.                 -11518334563.0 /  33831184612500.0,
  301.                 1912306948.0 /  13532473845.0
  302.             }, {
  303.                 -13613986967.0 / 31741908048.0,
  304.                 0.0,
  305.                 0.0,
  306.                 0.0,
  307.                 0.0,
  308.                 -4755612631.0 / 1012344804.0,
  309.                 42939257944576.0 / 5588559685701.0,
  310.                 77881972900277.0 / 19140370552944.0,
  311.                 22719829234375.0 / 63689648654052.0,
  312.                 0.0,
  313.                 0.0,
  314.                 0.0,
  315.                 -1199007803.0 / 857031517296.0,
  316.                 157882067000.0 / 53564469831.0,
  317.                 -290468882375.0 / 31741908048.0
  318.             }
  319.         };
  320.     }

  322.     /** {@inheritDoc} */
  323.     @Override
  324.     public double[] getB() {
  325.         return new double[] {
  326.             104257.0/1920240.0,
  327.             0.0,
  328.             0.0,
  329.             0.0,
  330.             0.0,
  331.             3399327.0/763840.0,
  332.             66578432.0/35198415.0,
  333.             -1674902723.0/288716400.0,
  334.             54980371265625.0/176692375811392.0,
  335.             -734375.0/4826304.0,
  336.             171414593.0/851261400.0,
  337.             137909.0/3084480.0,
  338.             0.0,
  339.             0.0,
  340.             0.0,
  341.             0.0
  342.         };
  343.     }

  345.     /** {@inheritDoc} */
  346.     @Override
  347.     protected DormandPrince853StateInterpolator
  348.     createInterpolator(final boolean forward, double[][] yDotK,
  349.                        final ODEStateAndDerivative globalPreviousState,
  350.                        final ODEStateAndDerivative globalCurrentState,
  351.                        final EquationsMapper mapper) {
  352.         return new DormandPrince853StateInterpolator(forward, yDotK,
  353.                                                     globalPreviousState, globalCurrentState,
  354.                                                     globalPreviousState, globalCurrentState,
  355.                                                     mapper);
  356.     }

  358.     /** {@inheritDoc} */
  359.     @Override
  360.     public int getOrder() {
  361.         return 8;
  362.     }

  364.     /** {@inheritDoc} */
  365.     @Override
  366.     protected double estimateError(final double[][] yDotK,
  367.                                    final double[] y0, final double[] y1,
  368.                                    final double h) {

  370.         final StepsizeHelper helper = getStepSizeHelper();
  371.         double error1 = 0;
  372.         double error2 = 0;

  374.         for (int j = 0; j < helper.getMainSetDimension(); ++j) {
  375.             final double errSum1 = E1_01 * yDotK[0][j]  + E1_06 * yDotK[5][j] +
  376.                                    E1_07 * yDotK[6][j]  + E1_08 * yDotK[7][j] +
  377.                                    E1_09 * yDotK[8][j]  + E1_10 * yDotK[9][j] +
  378.                                    E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
  379.             final double errSum2 = E2_01 * yDotK[0][j]  + E2_06 * yDotK[5][j] +
  380.                                    E2_07 * yDotK[6][j]  + E2_08 * yDotK[7][j] +
  381.                                    E2_09 * yDotK[8][j]  + E2_10 * yDotK[9][j] +
  382.                                    E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];

  384.             final double tol = helper.getTolerance(j, FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j])));
  385.             final double ratio1  = errSum1 / tol;
  386.             error1        += ratio1 * ratio1;
  387.             final double ratio2  = errSum2 / tol;
  388.             error2        += ratio2 * ratio2;
  389.         }

  391.         double den = error1 + 0.01 * error2;
  392.         if (den <= 0.0) {
  393.             den = 1.0;
  394.         }

  396.         return FastMath.abs(h) * error1 / FastMath.sqrt(helper.getMainSetDimension() * den);

  398.     }

  400. }
Created with JaCoCo 0.8.12.202403310830