DormandPrince853FieldIntegrator.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.legacy.core.Field;
  21. import org.apache.commons.math4.legacy.core.RealFieldElement;
  22. import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
  23. import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
  24. import org.apache.commons.math4.legacy.core.MathArrays;


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

  61. public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
  62.     extends EmbeddedRungeKuttaFieldIntegrator<T> {

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

  67.     /** First error weights array, element 1. */
  68.     private final T e1_01;

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

  72.     /** First error weights array, element 6. */
  73.     private final T e1_06;

  75.     /** First error weights array, element 7. */
  76.     private final T e1_07;

  78.     /** First error weights array, element 8. */
  79.     private final T e1_08;

  81.     /** First error weights array, element 9. */
  82.     private final T e1_09;

  84.     /** First error weights array, element 10. */
  85.     private final T e1_10;

  87.     /** First error weights array, element 11. */
  88.     private final T e1_11;

  90.     /** First error weights array, element 12. */
  91.     private final T e1_12;


  94.     /** Second error weights array, element 1. */
  95.     private final T e2_01;

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

  99.     /** Second error weights array, element 6. */
  100.     private final T e2_06;

  102.     /** Second error weights array, element 7. */
  103.     private final T e2_07;

  105.     /** Second error weights array, element 8. */
  106.     private final T e2_08;

  108.     /** Second error weights array, element 9. */
  109.     private final T e2_09;

  111.     /** Second error weights array, element 10. */
  112.     private final T e2_10;

  114.     /** Second error weights array, element 11. */
  115.     private final T e2_11;

  117.     /** Second error weights array, element 12. */
  118.     private final T e2_12;

  120.     /** Simple constructor.
  121.      * Build an eighth order Dormand-Prince integrator with the given step bounds
  122.      * @param field field to which the time and state vector elements belong
  123.      * @param minStep minimal step (sign is irrelevant, regardless of
  124.      * integration direction, forward or backward), the last step can
  125.      * be smaller than this
  126.      * @param maxStep maximal step (sign is irrelevant, regardless of
  127.      * integration direction, forward or backward), the last step can
  128.      * be smaller than this
  129.      * @param scalAbsoluteTolerance allowed absolute error
  130.      * @param scalRelativeTolerance allowed relative error
  131.      */
  132.     public DormandPrince853FieldIntegrator(final Field<T> field,
  133.                                            final double minStep, final double maxStep,
  134.                                            final double scalAbsoluteTolerance,
  135.                                            final double scalRelativeTolerance) {
  136.         super(field, METHOD_NAME, 12,
  137.               minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
  138.         e1_01 = fraction(        116092271.0,       8848465920.0);
  139.         e1_06 = fraction(         -1871647.0,          1527680.0);
  140.         e1_07 = fraction(        -69799717.0,        140793660.0);
  141.         e1_08 = fraction(    1230164450203.0,     739113984000.0);
  142.         e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
  143.         e1_10 = fraction(        464500805.0,       1389975552.0);
  144.         e1_11 = fraction(    1606764981773.0,   19613062656000.0);
  145.         e1_12 = fraction(          -137909.0,          6168960.0);
  146.         e2_01 = fraction(          -364463.0,          1920240.0);
  147.         e2_06 = fraction(          3399327.0,           763840.0);
  148.         e2_07 = fraction(         66578432.0,         35198415.0);
  149.         e2_08 = fraction(      -1674902723.0,        288716400.0);
  150.         e2_09 = fraction(  -74684743568175.0,  176692375811392.0);
  151.         e2_10 = fraction(          -734375.0,          4826304.0);
  152.         e2_11 = fraction(        171414593.0,        851261400.0);
  153.         e2_12 = fraction(            69869.0,          3084480.0);
  154.     }

  156.     /** Simple constructor.
  157.      * Build an eighth order Dormand-Prince integrator with the given step bounds
  158.      * @param field field to which the time and state vector elements belong
  159.      * @param minStep minimal step (sign is irrelevant, regardless of
  160.      * integration direction, forward or backward), the last step can
  161.      * be smaller than this
  162.      * @param maxStep maximal step (sign is irrelevant, regardless of
  163.      * integration direction, forward or backward), the last step can
  164.      * be smaller than this
  165.      * @param vecAbsoluteTolerance allowed absolute error
  166.      * @param vecRelativeTolerance allowed relative error
  167.      */
  168.     public DormandPrince853FieldIntegrator(final Field<T> field,
  169.                                            final double minStep, final double maxStep,
  170.                                            final double[] vecAbsoluteTolerance,
  171.                                            final double[] vecRelativeTolerance) {
  172.         super(field, METHOD_NAME, 12,
  173.               minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  174.         e1_01 = fraction(        116092271.0,       8848465920.0);
  175.         e1_06 = fraction(         -1871647.0,          1527680.0);
  176.         e1_07 = fraction(        -69799717.0,        140793660.0);
  177.         e1_08 = fraction(    1230164450203.0,     739113984000.0);
  178.         e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
  179.         e1_10 = fraction(        464500805.0,       1389975552.0);
  180.         e1_11 = fraction(    1606764981773.0,   19613062656000.0);
  181.         e1_12 = fraction(          -137909.0,          6168960.0);
  182.         e2_01 = fraction(          -364463.0,          1920240.0);
  183.         e2_06 = fraction(          3399327.0,           763840.0);
  184.         e2_07 = fraction(         66578432.0,         35198415.0);
  185.         e2_08 = fraction(      -1674902723.0,        288716400.0);
  186.         e2_09 = fraction(  -74684743568175.0,  176692375811392.0);
  187.         e2_10 = fraction(          -734375.0,          4826304.0);
  188.         e2_11 = fraction(        171414593.0,        851261400.0);
  189.         e2_12 = fraction(            69869.0,          3084480.0);
  190.     }

  192.     /** {@inheritDoc} */
  193.     @Override
  194.     public T[] getC() {

  196.         final T sqrt6 = getField().getOne().multiply(6).sqrt();

  198.         final T[] c = MathArrays.buildArray(getField(), 15);
  199.         c[ 0] = sqrt6.add(-6).divide(-67.5);
  200.         c[ 1] = sqrt6.add(-6).divide(-45.0);
  201.         c[ 2] = sqrt6.add(-6).divide(-30.0);
  202.         c[ 3] = sqrt6.add( 6).divide( 30.0);
  203.         c[ 4] = fraction(1, 3);
  204.         c[ 5] = fraction(1, 4);
  205.         c[ 6] = fraction(4, 13);
  206.         c[ 7] = fraction(127, 195);
  207.         c[ 8] = fraction(3, 5);
  208.         c[ 9] = fraction(6, 7);
  209.         c[10] = getField().getOne();
  210.         c[11] = getField().getOne();
  211.         c[12] = fraction(1.0, 10.0);
  212.         c[13] = fraction(1.0, 5.0);
  213.         c[14] = fraction(7.0, 9.0);

  215.         return c;
  216.     }

  218.     /** {@inheritDoc} */
  219.     @Override
  220.     public T[][] getA() {

  222.         final T sqrt6 = getField().getOne().multiply(6).sqrt();

  224.         final T[][] a = MathArrays.buildArray(getField(), 15, -1);
  225.         for (int i = 0; i < a.length; ++i) {
  226.             a[i] = MathArrays.buildArray(getField(), i + 1);
  227.         }

  229.         a[ 0][ 0] = sqrt6.add(-6).divide(-67.5);

  231.         a[ 1][ 0] = sqrt6.add(-6).divide(-180);
  232.         a[ 1][ 1] = sqrt6.add(-6).divide( -60);

  234.         a[ 2][ 0] = sqrt6.add(-6).divide(-120);
  235.         a[ 2][ 1] = getField().getZero();
  236.         a[ 2][ 2] = sqrt6.add(-6).divide( -40);

  238.         a[ 3][ 0] = sqrt6.multiply(107).add(462).divide( 3000);
  239.         a[ 3][ 1] = getField().getZero();
  240.         a[ 3][ 2] = sqrt6.multiply(197).add(402).divide(-1000);
  241.         a[ 3][ 3] = sqrt6.multiply( 73).add(168).divide(  375);

  243.         a[ 4][ 0] = fraction(1, 27);
  244.         a[ 4][ 1] = getField().getZero();
  245.         a[ 4][ 2] = getField().getZero();
  246.         a[ 4][ 3] = sqrt6.add( 16).divide( 108);
  247.         a[ 4][ 4] = sqrt6.add(-16).divide(-108);

  249.         a[ 5][ 0] = fraction(19, 512);
  250.         a[ 5][ 1] = getField().getZero();
  251.         a[ 5][ 2] = getField().getZero();
  252.         a[ 5][ 3] = sqrt6.multiply( 23).add(118).divide(1024);
  253.         a[ 5][ 4] = sqrt6.multiply(-23).add(118).divide(1024);
  254.         a[ 5][ 5] = fraction(-9, 512);

  256.         a[ 6][ 0] = fraction(13772, 371293);
  257.         a[ 6][ 1] = getField().getZero();
  258.         a[ 6][ 2] = getField().getZero();
  259.         a[ 6][ 3] = sqrt6.multiply( 4784).add(51544).divide(371293);
  260.         a[ 6][ 4] = sqrt6.multiply(-4784).add(51544).divide(371293);
  261.         a[ 6][ 5] = fraction(-5688, 371293);
  262.         a[ 6][ 6] = fraction( 3072, 371293);

  264.         a[ 7][ 0] = fraction(58656157643.0, 93983540625.0);
  265.         a[ 7][ 1] = getField().getZero();
  266.         a[ 7][ 2] = getField().getZero();
  267.         a[ 7][ 3] = sqrt6.multiply(-318801444819.0).add(-1324889724104.0).divide(626556937500.0);
  268.         a[ 7][ 4] = sqrt6.multiply( 318801444819.0).add(-1324889724104.0).divide(626556937500.0);
  269.         a[ 7][ 5] = fraction(96044563816.0, 3480871875.0);
  270.         a[ 7][ 6] = fraction(5682451879168.0, 281950621875.0);
  271.         a[ 7][ 7] = fraction(-165125654.0, 3796875.0);

  273.         a[ 8][ 0] = fraction(8909899.0, 18653125.0);
  274.         a[ 8][ 1] = getField().getZero();
  275.         a[ 8][ 2] = getField().getZero();
  276.         a[ 8][ 3] = sqrt6.multiply(-1137963.0).add(-4521408.0).divide(2937500.0);
  277.         a[ 8][ 4] = sqrt6.multiply( 1137963.0).add(-4521408.0).divide(2937500.0);
  278.         a[ 8][ 5] = fraction(96663078.0, 4553125.0);
  279.         a[ 8][ 6] = fraction(2107245056.0, 137915625.0);
  280.         a[ 8][ 7] = fraction(-4913652016.0, 147609375.0);
  281.         a[ 8][ 8] = fraction(-78894270.0, 3880452869.0);

  283.         a[ 9][ 0] = fraction(-20401265806.0, 21769653311.0);
  284.         a[ 9][ 1] = getField().getZero();
  285.         a[ 9][ 2] = getField().getZero();
  286.         a[ 9][ 3] = sqrt6.multiply( 94326.0).add(354216.0).divide(112847.0);
  287.         a[ 9][ 4] = sqrt6.multiply(-94326.0).add(354216.0).divide(112847.0);
  288.         a[ 9][ 5] = fraction(-43306765128.0, 5313852383.0);
  289.         a[ 9][ 6] = fraction(-20866708358144.0, 1126708119789.0);
  290.         a[ 9][ 7] = fraction(14886003438020.0, 654632330667.0);
  291.         a[ 9][ 8] = fraction(35290686222309375.0, 14152473387134411.0);
  292.         a[ 9][ 9] = fraction(-1477884375.0, 485066827.0);

  294.         a[10][ 0] = fraction(39815761.0, 17514443.0);
  295.         a[10][ 1] = getField().getZero();
  296.         a[10][ 2] = getField().getZero();
  297.         a[10][ 3] = sqrt6.multiply(-960905.0).add(-3457480.0).divide(551636.0);
  298.         a[10][ 4] = sqrt6.multiply( 960905.0).add(-3457480.0).divide(551636.0);
  299.         a[10][ 5] = fraction(-844554132.0, 47026969.0);
  300.         a[10][ 6] = fraction(8444996352.0, 302158619.0);
  301.         a[10][ 7] = fraction(-2509602342.0, 877790785.0);
  302.         a[10][ 8] = fraction(-28388795297996250.0, 3199510091356783.0);
  303.         a[10][ 9] = fraction(226716250.0, 18341897.0);
  304.         a[10][10] = fraction(1371316744.0, 2131383595.0);

  306.         // the following stage is both for interpolation and the first stage in next step
  307.         // (the coefficients are identical to the B array)
  308.         a[11][ 0] = fraction(104257.0, 1920240.0);
  309.         a[11][ 1] = getField().getZero();
  310.         a[11][ 2] = getField().getZero();
  311.         a[11][ 3] = getField().getZero();
  312.         a[11][ 4] = getField().getZero();
  313.         a[11][ 5] = fraction(3399327.0, 763840.0);
  314.         a[11][ 6] = fraction(66578432.0, 35198415.0);
  315.         a[11][ 7] = fraction(-1674902723.0, 288716400.0);
  316.         a[11][ 8] = fraction(54980371265625.0, 176692375811392.0);
  317.         a[11][ 9] = fraction(-734375.0, 4826304.0);
  318.         a[11][10] = fraction(171414593.0, 851261400.0);
  319.         a[11][11] = fraction(137909.0, 3084480.0);

  321.         // the following stages are for interpolation only
  322.         a[12][ 0] = fraction(      13481885573.0, 240030000000.0);
  323.         a[12][ 1] = getField().getZero();
  324.         a[12][ 2] = getField().getZero();
  325.         a[12][ 3] = getField().getZero();
  326.         a[12][ 4] = getField().getZero();
  327.         a[12][ 5] = getField().getZero();
  328.         a[12][ 6] = fraction(     139418837528.0, 549975234375.0);
  329.         a[12][ 7] = fraction(  -11108320068443.0, 45111937500000.0);
  330.         a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0);
  331.         a[12][ 9] = fraction(         57799439.0, 377055000.0);
  332.         a[12][10] = fraction(     793322643029.0, 96734250000000.0);
  333.         a[12][11] = fraction(       1458939311.0, 192780000000.0);
  334.         a[12][12]  = fraction(            -4149.0, 500000.0);

  336.         a[13][ 0] = fraction(    1595561272731.0, 50120273500000.0);
  337.         a[13][ 1] = getField().getZero();
  338.         a[13][ 2] = getField().getZero();
  339.         a[13][ 3] = getField().getZero();
  340.         a[13][ 4] = getField().getZero();
  341.         a[13][ 5] = fraction(     975183916491.0, 34457688031250.0);
  342.         a[13][ 6] = fraction(   38492013932672.0, 718912673015625.0);
  343.         a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0);
  344.         a[13][ 8] = getField().getZero();
  345.         a[13][ 9] = getField().getZero();
  346.         a[13][10] = fraction(   -2538710946863.0, 23431227861250000.0);
  347.         a[13][11] = fraction(       8824659001.0, 23066716781250.0);
  348.         a[13][12] = fraction(     -11518334563.0, 33831184612500.0);
  349.         a[13][13] = fraction(       1912306948.0, 13532473845.0);

  351.         a[14][ 0] = fraction(     -13613986967.0, 31741908048.0);
  352.         a[14][ 1] = getField().getZero();
  353.         a[14][ 2] = getField().getZero();
  354.         a[14][ 3] = getField().getZero();
  355.         a[14][ 4] = getField().getZero();
  356.         a[14][ 5] = fraction(      -4755612631.0, 1012344804.0);
  357.         a[14][ 6] = fraction(   42939257944576.0, 5588559685701.0);
  358.         a[14][ 7] = fraction(   77881972900277.0, 19140370552944.0);
  359.         a[14][ 8] = fraction(   22719829234375.0, 63689648654052.0);
  360.         a[14][ 9] = getField().getZero();
  361.         a[14][10] = getField().getZero();
  362.         a[14][11] = getField().getZero();
  363.         a[14][12] = fraction(      -1199007803.0, 857031517296.0);
  364.         a[14][13] = fraction(     157882067000.0, 53564469831.0);
  365.         a[14][14] = fraction(    -290468882375.0, 31741908048.0);

  367.         return a;
  368.     }

  370.     /** {@inheritDoc} */
  371.     @Override
  372.     public T[] getB() {
  373.         final T[] b = MathArrays.buildArray(getField(), 16);
  374.         b[ 0] = fraction(104257, 1920240);
  375.         b[ 1] = getField().getZero();
  376.         b[ 2] = getField().getZero();
  377.         b[ 3] = getField().getZero();
  378.         b[ 4] = getField().getZero();
  379.         b[ 5] = fraction(        3399327.0,          763840.0);
  380.         b[ 6] = fraction(       66578432.0,        35198415.0);
  381.         b[ 7] = fraction(    -1674902723.0,       288716400.0);
  382.         b[ 8] = fraction( 54980371265625.0, 176692375811392.0);
  383.         b[ 9] = fraction(        -734375.0,         4826304.0);
  384.         b[10] = fraction(      171414593.0,       851261400.0);
  385.         b[11] = fraction(         137909.0,         3084480.0);
  386.         b[12] = getField().getZero();
  387.         b[13] = getField().getZero();
  388.         b[14] = getField().getZero();
  389.         b[15] = getField().getZero();
  390.         return b;
  391.     }

  393.     /** {@inheritDoc} */
  394.     @Override
  395.     protected DormandPrince853FieldStepInterpolator<T>
  396.         createInterpolator(final boolean forward, T[][] yDotK,
  397.                            final FieldODEStateAndDerivative<T> globalPreviousState,
  398.                            final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
  399.         return new DormandPrince853FieldStepInterpolator<>(getField(), forward, yDotK,
  400.                                                             globalPreviousState, globalCurrentState,
  401.                                                             globalPreviousState, globalCurrentState,
  402.                                                             mapper);
  403.     }

  405.     /** {@inheritDoc} */
  406.     @Override
  407.     public int getOrder() {
  408.         return 8;
  409.     }

  411.     /** {@inheritDoc} */
  412.     @Override
  413.     protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
  414.         T error1 = h.getField().getZero();
  415.         T error2 = h.getField().getZero();

  417.         for (int j = 0; j < mainSetDimension; ++j) {
  418.             final T errSum1 =      yDotK[ 0][j].multiply(e1_01).
  419.                                add(yDotK[ 5][j].multiply(e1_06)).
  420.                                add(yDotK[ 6][j].multiply(e1_07)).
  421.                                add(yDotK[ 7][j].multiply(e1_08)).
  422.                                add(yDotK[ 8][j].multiply(e1_09)).
  423.                                add(yDotK[ 9][j].multiply(e1_10)).
  424.                                add(yDotK[10][j].multiply(e1_11)).
  425.                                add(yDotK[11][j].multiply(e1_12));
  426.             final T errSum2 =      yDotK[ 0][j].multiply(e2_01).
  427.                                add(yDotK[ 5][j].multiply(e2_06)).
  428.                                add(yDotK[ 6][j].multiply(e2_07)).
  429.                                add(yDotK[ 7][j].multiply(e2_08)).
  430.                                add(yDotK[ 8][j].multiply(e2_09)).
  431.                                add(yDotK[ 9][j].multiply(e2_10)).
  432.                                add(yDotK[10][j].multiply(e2_11)).
  433.                                add(yDotK[11][j].multiply(e2_12));

  435.             final T yScale = RealFieldElement.max(y0[j].abs(), y1[j].abs());
  436.             final T tol = vecAbsoluteTolerance == null ?
  437.                           yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
  438.                           yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]);
  439.             final T ratio1  = errSum1.divide(tol);
  440.             error1        = error1.add(ratio1.multiply(ratio1));
  441.             final T ratio2  = errSum2.divide(tol);
  442.             error2        = error2.add(ratio2.multiply(ratio2));
  443.         }

  445.         T den = error1.add(error2.multiply(0.01));
  446.         if (den.getReal() <= 0.0) {
  447.             den = h.getField().getOne();
  448.         }

  450.         return h.abs().multiply(error1).divide(den.multiply(mainSetDimension).sqrt());
  451.     }
  452. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT