DormandPrince853StateInterpolator.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.interpolators;

  20. import org.hipparchus.ode.EquationsMapper;
  21. import org.hipparchus.ode.ODEStateAndDerivative;
  22. import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;

  24. /**
  25.  * This class represents an interpolator over the last step during an
  26.  * ODE integration for the 8(5,3) Dormand-Prince integrator.
  27.  *
  28.  * @see DormandPrince853Integrator
  29.  *
  30.  */

  32. public class DormandPrince853StateInterpolator extends RungeKuttaStateInterpolator {

  34.     /** Serializable version identifier. */
  35.     private static final long serialVersionUID = 20160328L;

  37.     /** Interpolation weights. */
  38.     private static final double[][] D = {
  39.         {
  40.             // this row is the same as the b array
  41.                     104257.0 / 1920240.0,
  42.                          0.0,
  43.                          0.0,
  44.                          0.0,
  45.                          0.0,
  46.                    3399327.0 / 763840.0,
  47.                   66578432.0 / 35198415.0,
  48.                -1674902723.0 / 288716400.0,
  49.             54980371265625.0 / 176692375811392.0,
  50.                    -734375.0 / 4826304.0,
  51.                  171414593.0 / 851261400.0,
  52.                     137909.0 / 3084480.0,
  53.                          0.0,
  54.                          0.0,
  55.                          0.0,
  56.                          0.0,
  57.         }, {
  58.                     1815983.0 / 1920240.0,
  59.                           0.0,
  60.                           0.0,
  61.                           0.0,
  62.                           0.0,
  63.                    -3399327.0 / 763840.0,
  64.                   -66578432.0 / 35198415.0,
  65.                  1674902723.0 / 288716400.0,
  66.             -54980371265625.0 / 176692375811392.0,
  67.                      734375.0 / 4826304.0,
  68.                  -171414593.0 / 851261400.0,
  69.                     -137909.0 / 3084480.0,
  70.                           0.0,
  71.                           0.0,
  72.                           0.0,
  73.                           0.0,
  74.         }, {
  75.                     -855863.0 / 960120.0,
  76.                           0.0,
  77.                           0.0,
  78.                           0.0,
  79.                           0.0,
  80.                     3399327.0 / 381920.0,
  81.                   133156864.0 / 35198415.0,
  82.                 -1674902723.0 / 144358200.0,
  83.              54980371265625.0 / 88346187905696.0,
  84.                     -734375.0 / 2413152.0,
  85.                   171414593.0 / 425630700.0,
  86.                      137909.0 / 1542240.0,
  87.                          -1.0,
  88.                           0.0,
  89.                           0.0,
  90.                           0.0
  91.         }, {
  92.                -17751989329.0 / 2106076560.0,
  93.                           0.0,
  94.                           0.0,
  95.                           0.0,
  96.                           0.0,
  97.                  4272954039.0 / 7539864640.0,
  98.               -118476319744.0 / 38604839385.0,
  99.                755123450731.0 / 316657731600.0,
  100.         3692384461234828125.0 / 1744130441634250432.0,
  101.                 -4612609375.0 / 5293382976.0,
  102.               2091772278379.0 / 933644586600.0,
  103.                  2136624137.0 / 3382989120.0,
  104.                     -126493.0 / 1421424.0,
  105.                    98350000.0 / 5419179.0,
  106.                   -18878125.0 / 2053168.0,
  107.                 -1944542619.0 / 438351368.0
  108.         }, {
  109.                 32941697297.0 / 3159114840.0,
  110.                           0.0,
  111.                           0.0,
  112.                           0.0,
  113.                           0.0,
  114.                456696183123.0 / 1884966160.0,
  115.              19132610714624.0 / 115814518155.0,
  116.            -177904688592943.0 / 474986597400.0,
  117.        -4821139941836765625.0 / 218016305204281304.0,
  118.                 30702015625.0 / 3970037232.0,
  119.             -85916079474274.0 / 2800933759800.0,
  120.                 -5919468007.0 / 634310460.0,
  121.                     2479159.0 / 157936.0,
  122.                   -18750000.0 / 602131.0,
  123.                   -19203125.0 / 2053168.0,
  124.                 15700361463.0 / 438351368.0
  125.         }, {
  126.                 12627015655.0 / 631822968.0,
  127.                           0.0,
  128.                           0.0,
  129.                           0.0,
  130.                           0.0,
  131.                -72955222965.0 / 188496616.0,
  132.             -13145744952320.0 / 69488710893.0,
  133.              30084216194513.0 / 56998391688.0,
  134.         -296858761006640625.0 / 25648977082856624.0,
  135.                   569140625.0 / 82709109.0,
  136.                -18684190637.0 / 18672891732.0,
  137.                    69644045.0 / 89549712.0,
  138.                   -11847025.0 / 4264272.0,
  139.                  -978650000.0 / 16257537.0,
  140.                   519371875.0 / 6159504.0,
  141.                  5256837225.0 / 438351368.0
  142.         }, {
  143.                  -450944925.0 / 17550638.0,
  144.                           0.0,
  145.                           0.0,
  146.                           0.0,
  147.                           0.0,
  148.                -14532122925.0 / 94248308.0,
  149.               -595876966400.0 / 2573655959.0,
  150.                188748653015.0 / 527762886.0,
  151.         2545485458115234375.0 / 27252038150535163.0,
  152.                 -1376953125.0 / 36759604.0,
  153.                 53995596795.0 / 518691437.0,
  154.                   210311225.0 / 7047894.0,
  155.                    -1718875.0 / 39484.0,
  156.                    58000000.0 / 602131.0,
  157.                    -1546875.0 / 39484.0,
  158.                 -1262172375.0 / 8429834.0
  159.         }
  160.     };

  162.     /** Simple constructor.
  163.      * @param forward integration direction indicator
  164.      * @param yDotK slopes at the intermediate points
  165.      * @param globalPreviousState start of the global step
  166.      * @param globalCurrentState end of the global step
  167.      * @param softPreviousState start of the restricted step
  168.      * @param softCurrentState end of the restricted step
  169.      * @param mapper equations mapper for the all equations
  170.      */
  171.     public DormandPrince853StateInterpolator(final boolean forward,
  172.                                              final double[][] yDotK,
  173.                                              final ODEStateAndDerivative globalPreviousState,
  174.                                              final ODEStateAndDerivative globalCurrentState,
  175.                                              final ODEStateAndDerivative softPreviousState,
  176.                                              final ODEStateAndDerivative softCurrentState,
  177.                                              final EquationsMapper mapper) {
  178.         super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
  179.     }

  181.     /** {@inheritDoc} */
  182.     @Override
  183.     protected DormandPrince853StateInterpolator create(final boolean newForward, final double[][] newYDotK,
  184.                                                        final ODEStateAndDerivative newGlobalPreviousState,
  185.                                                        final ODEStateAndDerivative newGlobalCurrentState,
  186.                                                        final ODEStateAndDerivative newSoftPreviousState,
  187.                                                        final ODEStateAndDerivative newSoftCurrentState,
  188.                                                        final EquationsMapper newMapper) {
  189.         return new DormandPrince853StateInterpolator(newForward, newYDotK,
  190.                                                      newGlobalPreviousState, newGlobalCurrentState,
  191.                                                      newSoftPreviousState, newSoftCurrentState,
  192.                                                      newMapper);
  193.     }

  195.     /** {@inheritDoc} */
  196.     @Override
  197.     protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
  198.                                                                            final double time, final double theta,
  199.                                                                            final double thetaH, final double oneMinusThetaH) {

  201.         final double eta      = 1.0 - theta;
  202.         final double twoTheta = 2 * theta;
  203.         final double theta2   = theta * theta;
  204.         final double dot1     = 1.0 - twoTheta;
  205.         final double dot2     = theta * (2 - 3 * theta);
  206.         final double dot3     = twoTheta * (theta * (twoTheta - 3) + 1);
  207.         final double dot4     = theta2   * (theta * (5 * theta - 8) + 3);
  208.         final double dot5     = theta2   * (theta * (theta * (15 - 6 * theta) - 12) + 3);
  209.         final double dot6     = theta2   * (theta * (theta * (theta * (18 - 7 * theta) - 15) + 4));
  210.         final double[] interpolatedState;
  211.         final double[] interpolatedDerivatives;


  214.         if (getGlobalPreviousState() != null && theta <= 0.5) {
  215.             final double f0 = thetaH;
  216.             final double f1 = f0 * eta;
  217.             final double f2 = f1 * theta;
  218.             final double f3 = f2 * eta;
  219.             final double f4 = f3 * theta;
  220.             final double f5 = f4 * eta;
  221.             final double f6 = f5 * theta;
  222.             final double[] p = new double[16];
  223.             final double[] q = new double[16];
  224.             for (int i = 0; i < p.length; ++i) {
  225.                 p[i] = f0 * D[0][i] +   f1 * D[1][i] +   f2 * D[2][i] +   f3 * D[3][i] +
  226.                        f4 * D[4][i] +   f5 * D[5][i] +   f6 * D[6][i];
  227.                 q[i] = D[0][i] + dot1 * D[1][i] + dot2 * D[2][i] + dot3 * D[3][i] +
  228.                                  dot4 * D[4][i] + dot5 * D[5][i] + dot6 * D[6][i];
  229.             }
  230.             interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
  231.                                                                      p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
  232.             interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
  233.                                                                   q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
  234.         } else {
  235.             final double f0 = -oneMinusThetaH;
  236.             final double f1 = -f0 * theta;
  237.             final double f2 =  f1 * theta;
  238.             final double f3 =  f2 * eta;
  239.             final double f4 =  f3 * theta;
  240.             final double f5 =  f4 * eta;
  241.             final double f6 =  f5 * theta;
  242.             final double[] p = new double[16];
  243.             final double[] q = new double[16];
  244.             for (int i = 0; i < p.length; ++i) {
  245.                 p[i] = f0 * D[0][i] +   f1 * D[1][i] +   f2 * D[2][i] +   f3 * D[3][i] +
  246.                        f4 * D[4][i] +   f5 * D[5][i] +   f6 * D[6][i];
  247.                 q[i] = D[0][i] + dot1 * D[1][i] + dot2 * D[2][i] + dot3 * D[3][i] +
  248.                                  dot4 * D[4][i] + dot5 * D[5][i] + dot6 * D[6][i];
  249.             }
  250.             interpolatedState       = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
  251.                                                                     p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
  252.             interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
  253.                                                                   q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
  254.         }

  256.         return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);

  258.     }

  260. }
Created with JaCoCo 0.8.12.202403310830