DormandPrince853FieldStateInterpolator.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.  *      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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.ode.nonstiff.interpolators;

  25. import org.hipparchus.CalculusFieldElement;
  26. import org.hipparchus.Field;
  27. import org.hipparchus.exception.MathIllegalStateException;
  28. import org.hipparchus.ode.FieldEquationsMapper;
  29. import org.hipparchus.ode.FieldODEStateAndDerivative;
  30. import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
  31. import org.hipparchus.util.MathArrays;

  33. /**
  34.  * This class represents an interpolator over the last step during an
  35.  * ODE integration for the 8(5,3) Dormand-Prince integrator.
  36.  *
  37.  * @see DormandPrince853FieldIntegrator
  38.  *
  39.  * @param <T> the type of the field elements
  40.  */

  42. public class DormandPrince853FieldStateInterpolator<T extends CalculusFieldElement<T>>
  43.     extends RungeKuttaFieldStateInterpolator<T> {

  45.     /** Interpolation weights.
  46.      * (beware that only the non-null values are in the table)
  47.      */
  48.     private final T[][] d;

  50.     /** Simple constructor.
  51.      * @param field field to which the time and state vector elements belong
  52.      * @param forward integration direction indicator
  53.      * @param yDotK slopes at the intermediate points
  54.      * @param globalPreviousState start of the global step
  55.      * @param globalCurrentState end of the global step
  56.      * @param softPreviousState start of the restricted step
  57.      * @param softCurrentState end of the restricted step
  58.      * @param mapper equations mapper for the all equations
  59.      */
  60.     public DormandPrince853FieldStateInterpolator(final Field<T> field, final boolean forward,
  61.                                                   final T[][] yDotK,
  62.                                                   final FieldODEStateAndDerivative<T> globalPreviousState,
  63.                                                   final FieldODEStateAndDerivative<T> globalCurrentState,
  64.                                                   final FieldODEStateAndDerivative<T> softPreviousState,
  65.                                                   final FieldODEStateAndDerivative<T> softCurrentState,
  66.                                                   final FieldEquationsMapper<T> mapper) {
  67.         super(field, forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
  68.                 mapper);
  69.         // interpolation weights
  70.         d = MathArrays.buildArray(field, 7, 16);

  72.         // this row is the same as the b array
  73.         d[0][ 0] = fraction(field, 104257, 1920240);
  74.         d[0][ 1] = field.getZero();
  75.         d[0][ 2] = field.getZero();
  76.         d[0][ 3] = field.getZero();
  77.         d[0][ 4] = field.getZero();
  78.         d[0][ 5] = fraction(field,         3399327.0,          763840.0);
  79.         d[0][ 6] = fraction(field,        66578432.0,        35198415.0);
  80.         d[0][ 7] = fraction(field,     -1674902723.0,       288716400.0);
  81.         d[0][ 8] = fraction(field,  54980371265625.0, 176692375811392.0);
  82.         d[0][ 9] = fraction(field,         -734375.0,         4826304.0);
  83.         d[0][10] = fraction(field,       171414593.0,       851261400.0);
  84.         d[0][11] = fraction(field,          137909.0,         3084480.0);
  85.         d[0][12] = field.getZero();
  86.         d[0][13] = field.getZero();
  87.         d[0][14] = field.getZero();
  88.         d[0][15] = field.getZero();

  90.         d[1][ 0] = d[0][ 0].negate().add(1);
  91.         d[1][ 1] = d[0][ 1].negate();
  92.         d[1][ 2] = d[0][ 2].negate();
  93.         d[1][ 3] = d[0][ 3].negate();
  94.         d[1][ 4] = d[0][ 4].negate();
  95.         d[1][ 5] = d[0][ 5].negate();
  96.         d[1][ 6] = d[0][ 6].negate();
  97.         d[1][ 7] = d[0][ 7].negate();
  98.         d[1][ 8] = d[0][ 8].negate();
  99.         d[1][ 9] = d[0][ 9].negate();
  100.         d[1][10] = d[0][10].negate();
  101.         d[1][11] = d[0][11].negate();
  102.         d[1][12] = d[0][12].negate(); // really 0
  103.         d[1][13] = d[0][13].negate(); // really 0
  104.         d[1][14] = d[0][14].negate(); // really 0
  105.         d[1][15] = d[0][15].negate(); // really 0

  107.         d[2][ 0] = d[0][ 0].multiply(2).subtract(1);
  108.         d[2][ 1] = d[0][ 1].multiply(2);
  109.         d[2][ 2] = d[0][ 2].multiply(2);
  110.         d[2][ 3] = d[0][ 3].multiply(2);
  111.         d[2][ 4] = d[0][ 4].multiply(2);
  112.         d[2][ 5] = d[0][ 5].multiply(2);
  113.         d[2][ 6] = d[0][ 6].multiply(2);
  114.         d[2][ 7] = d[0][ 7].multiply(2);
  115.         d[2][ 8] = d[0][ 8].multiply(2);
  116.         d[2][ 9] = d[0][ 9].multiply(2);
  117.         d[2][10] = d[0][10].multiply(2);
  118.         d[2][11] = d[0][11].multiply(2);
  119.         d[2][12] = d[0][12].multiply(2).subtract(1); // really -1
  120.         d[2][13] = d[0][13].multiply(2);             // really  0
  121.         d[2][14] = d[0][14].multiply(2);             // really  0
  122.         d[2][15] = d[0][15].multiply(2);             // really  0

  124.         d[3][ 0] = fraction(field,         -17751989329.0, 2106076560.0);
  125.         d[3][ 1] = field.getZero();
  126.         d[3][ 2] = field.getZero();
  127.         d[3][ 3] = field.getZero();
  128.         d[3][ 4] = field.getZero();
  129.         d[3][ 5] = fraction(field,           4272954039.0, 7539864640.0);
  130.         d[3][ 6] = fraction(field,        -118476319744.0, 38604839385.0);
  131.         d[3][ 7] = fraction(field,         755123450731.0, 316657731600.0);
  132.         d[3][ 8] = fraction(field,  3692384461234828125.0, 1744130441634250432.0);
  133.         d[3][ 9] = fraction(field,          -4612609375.0, 5293382976.0);
  134.         d[3][10] = fraction(field,        2091772278379.0, 933644586600.0);
  135.         d[3][11] = fraction(field,           2136624137.0, 3382989120.0);
  136.         d[3][12] = fraction(field,              -126493.0, 1421424.0);
  137.         d[3][13] = fraction(field,             98350000.0, 5419179.0);
  138.         d[3][14] = fraction(field,            -18878125.0, 2053168.0);
  139.         d[3][15] = fraction(field,          -1944542619.0, 438351368.0);

  141.         d[4][ 0] = fraction(field,          32941697297.0, 3159114840.0);
  142.         d[4][ 1] = field.getZero();
  143.         d[4][ 2] = field.getZero();
  144.         d[4][ 3] = field.getZero();
  145.         d[4][ 4] = field.getZero();
  146.         d[4][ 5] = fraction(field,         456696183123.0, 1884966160.0);
  147.         d[4][ 6] = fraction(field,       19132610714624.0, 115814518155.0);
  148.         d[4][ 7] = fraction(field,     -177904688592943.0, 474986597400.0);
  149.         d[4][ 8] = fraction(field, -4821139941836765625.0, 218016305204281304.0);
  150.         d[4][ 9] = fraction(field,          30702015625.0, 3970037232.0);
  151.         d[4][10] = fraction(field,      -85916079474274.0, 2800933759800.0);
  152.         d[4][11] = fraction(field,          -5919468007.0, 634310460.0);
  153.         d[4][12] = fraction(field,              2479159.0, 157936.0);
  154.         d[4][13] = fraction(field,            -18750000.0, 602131.0);
  155.         d[4][14] = fraction(field,            -19203125.0, 2053168.0);
  156.         d[4][15] = fraction(field,          15700361463.0, 438351368.0);

  158.         d[5][ 0] = fraction(field,          12627015655.0, 631822968.0);
  159.         d[5][ 1] = field.getZero();
  160.         d[5][ 2] = field.getZero();
  161.         d[5][ 3] = field.getZero();
  162.         d[5][ 4] = field.getZero();
  163.         d[5][ 5] = fraction(field,         -72955222965.0, 188496616.0);
  164.         d[5][ 6] = fraction(field,      -13145744952320.0, 69488710893.0);
  165.         d[5][ 7] = fraction(field,       30084216194513.0, 56998391688.0);
  166.         d[5][ 8] = fraction(field,  -296858761006640625.0, 25648977082856624.0);
  167.         d[5][ 9] = fraction(field,            569140625.0, 82709109.0);
  168.         d[5][10] = fraction(field,         -18684190637.0, 18672891732.0);
  169.         d[5][11] = fraction(field,             69644045.0, 89549712.0);
  170.         d[5][12] = fraction(field,            -11847025.0, 4264272.0);
  171.         d[5][13] = fraction(field,           -978650000.0, 16257537.0);
  172.         d[5][14] = fraction(field,            519371875.0, 6159504.0);
  173.         d[5][15] = fraction(field,           5256837225.0, 438351368.0);

  175.         d[6][ 0] = fraction(field,           -450944925.0, 17550638.0);
  176.         d[6][ 1] = field.getZero();
  177.         d[6][ 2] = field.getZero();
  178.         d[6][ 3] = field.getZero();
  179.         d[6][ 4] = field.getZero();
  180.         d[6][ 5] = fraction(field,         -14532122925.0, 94248308.0);
  181.         d[6][ 6] = fraction(field,        -595876966400.0, 2573655959.0);
  182.         d[6][ 7] = fraction(field,         188748653015.0, 527762886.0);
  183.         d[6][ 8] = fraction(field,  2545485458115234375.0, 27252038150535163.0);
  184.         d[6][ 9] = fraction(field,          -1376953125.0, 36759604.0);
  185.         d[6][10] = fraction(field,          53995596795.0, 518691437.0);
  186.         d[6][11] = fraction(field,            210311225.0, 7047894.0);
  187.         d[6][12] = fraction(field,             -1718875.0, 39484.0);
  188.         d[6][13] = fraction(field,             58000000.0, 602131.0);
  189.         d[6][14] = fraction(field,             -1546875.0, 39484.0);
  190.         d[6][15] = fraction(field,          -1262172375.0, 8429834.0);

  192.     }

  194.     /** {@inheritDoc} */
  195.     @Override
  196.     protected DormandPrince853FieldStateInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
  197.                                                                final FieldODEStateAndDerivative<T> newGlobalPreviousState,
  198.                                                                final FieldODEStateAndDerivative<T> newGlobalCurrentState,
  199.                                                                final FieldODEStateAndDerivative<T> newSoftPreviousState,
  200.                                                                final FieldODEStateAndDerivative<T> newSoftCurrentState,
  201.                                                                final FieldEquationsMapper<T> newMapper) {
  202.         return new DormandPrince853FieldStateInterpolator<>(newField, newForward, newYDotK,
  203.                                                              newGlobalPreviousState, newGlobalCurrentState,
  204.                                                              newSoftPreviousState, newSoftCurrentState,
  205.                                                              newMapper);
  206.     }

  208.     /** Create a fraction.
  209.      * @param field field to which the elements belong
  210.      * @param p numerator
  211.      * @param q denominator
  212.      * @return p/q computed in the instance field
  213.      */
  214.     private T fraction(final Field<T> field, final double p, final double q) {
  215.         return field.getZero().add(p).divide(q);
  216.     }

  218.     /** {@inheritDoc} */
  219.     @SuppressWarnings("unchecked")
  220.     @Override
  221.     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
  222.                                                                                    final T time, final T theta,
  223.                                                                                    final T thetaH, final T oneMinusThetaH)
  224.         throws MathIllegalStateException {

  226.         final T one      = time.getField().getOne();
  227.         final T eta      = one.subtract(theta);
  228.         final T twoTheta = theta.multiply(2);
  229.         final T theta2   = theta.multiply(theta);
  230.         final T dot1     = one.subtract(twoTheta);
  231.         final T dot2     = theta.multiply(theta.multiply(-3).add(2));
  232.         final T dot3     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
  233.         final T dot4     = theta2.multiply(theta.multiply(theta.multiply(5).subtract(8)).add(3));
  234.         final T dot5     = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(-6).add(15)).subtract(12)).add(3));
  235.         final T dot6     = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-7).add(18)).subtract(15)).add(4)));
  236.         final T[] interpolatedState;
  237.         final T[] interpolatedDerivatives;


  240.         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
  241.             final T f0 = thetaH;
  242.             final T f1 = f0.multiply(eta);
  243.             final T f2 = f1.multiply(theta);
  244.             final T f3 = f2.multiply(eta);
  245.             final T f4 = f3.multiply(theta);
  246.             final T f5 = f4.multiply(eta);
  247.             final T f6 = f5.multiply(theta);
  248.             final T[] p = MathArrays.buildArray(time.getField(), 16);
  249.             final T[] q = MathArrays.buildArray(time.getField(), 16);
  250.             for (int i = 0; i < p.length; ++i) {
  251.                 p[i] =     f0.multiply(d[0][i]).
  252.                        add(f1.multiply(d[1][i])).
  253.                        add(f2.multiply(d[2][i])).
  254.                        add(f3.multiply(d[3][i])).
  255.                        add(f4.multiply(d[4][i])).
  256.                        add(f5.multiply(d[5][i])).
  257.                        add(f6.multiply(d[6][i]));
  258.                 q[i] =                    d[0][i].
  259.                         add(dot1.multiply(d[1][i])).
  260.                         add(dot2.multiply(d[2][i])).
  261.                         add(dot3.multiply(d[3][i])).
  262.                         add(dot4.multiply(d[4][i])).
  263.                         add(dot5.multiply(d[5][i])).
  264.                         add(dot6.multiply(d[6][i]));
  265.             }
  266.             interpolatedState       = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
  267.                                                                      p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
  268.             interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
  269.                                                                   q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
  270.         } else {
  271.             final T f0 = oneMinusThetaH.negate();
  272.             final T f1 = f0.multiply(theta).negate();
  273.             final T f2 = f1.multiply(theta);
  274.             final T f3 = f2.multiply(eta);
  275.             final T f4 = f3.multiply(theta);
  276.             final T f5 = f4.multiply(eta);
  277.             final T f6 = f5.multiply(theta);
  278.             final T[] p = MathArrays.buildArray(time.getField(), 16);
  279.             final T[] q = MathArrays.buildArray(time.getField(), 16);
  280.             for (int i = 0; i < p.length; ++i) {
  281.                 p[i] =     f0.multiply(d[0][i]).
  282.                        add(f1.multiply(d[1][i])).
  283.                        add(f2.multiply(d[2][i])).
  284.                        add(f3.multiply(d[3][i])).
  285.                        add(f4.multiply(d[4][i])).
  286.                        add(f5.multiply(d[5][i])).
  287.                        add(f6.multiply(d[6][i]));
  288.                 q[i] =                    d[0][i].
  289.                         add(dot1.multiply(d[1][i])).
  290.                         add(dot2.multiply(d[2][i])).
  291.                         add(dot3.multiply(d[3][i])).
  292.                         add(dot4.multiply(d[4][i])).
  293.                         add(dot5.multiply(d[5][i])).
  294.                         add(dot6.multiply(d[6][i]));
  295.             }
  296.             interpolatedState       = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
  297.                                                                     p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
  298.             interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
  299.                                                                   q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
  300.         }

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

  304.     }

  306. }
Created with JaCoCo 0.8.12.202403310830