DormandPrince54StateInterpolator.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.DormandPrince54Integrator;

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

  32. public class DormandPrince54StateInterpolator extends RungeKuttaStateInterpolator {

  34.     /** Last row of the Butcher-array internal weights, element 0. */
  35.     private static final double A70 =    35.0 /  384.0;

  37.     // element 1 is zero, so it is neither stored nor used

  39.     /** Last row of the Butcher-array internal weights, element 2. */
  40.     private static final double A72 =   500.0 / 1113.0;

  42.     /** Last row of the Butcher-array internal weights, element 3. */
  43.     private static final double A73 =   125.0 /  192.0;

  45.     /** Last row of the Butcher-array internal weights, element 4. */
  46.     private static final double A74 = -2187.0 / 6784.0;

  48.     /** Last row of the Butcher-array internal weights, element 5. */
  49.     private static final double A75 =    11.0 /   84.0;

  51.     /** Shampine (1986) Dense output, element 0. */
  52.     private static final double D0 =  -12715105075.0 /  11282082432.0;

  54.     // element 1 is zero, so it is neither stored nor used

  56.     /** Shampine (1986) Dense output, element 2. */
  57.     private static final double D2 =   87487479700.0 /  32700410799.0;

  59.     /** Shampine (1986) Dense output, element 3. */
  60.     private static final double D3 =  -10690763975.0 /   1880347072.0;

  62.     /** Shampine (1986) Dense output, element 4. */
  63.     private static final double D4 =  701980252875.0 / 199316789632.0;

  65.     /** Shampine (1986) Dense output, element 5. */
  66.     private static final double D5 =   -1453857185.0 /    822651844.0;

  68.     /** Shampine (1986) Dense output, element 6. */
  69.     private static final double D6 =      69997945.0 /     29380423.0;

  71.     /** Serializable version identifier. */
  72.     private static final long serialVersionUID = 20160328L;

  74.     /** Simple constructor.
  75.      * @param forward integration direction indicator
  76.      * @param yDotK slopes at the intermediate points
  77.      * @param globalPreviousState start of the global step
  78.      * @param globalCurrentState end of the global step
  79.      * @param softPreviousState start of the restricted step
  80.      * @param softCurrentState end of the restricted step
  81.      * @param mapper equations mapper for the all equations
  82.      */
  83.     public DormandPrince54StateInterpolator(final boolean forward,
  84.                                      final double[][] yDotK,
  85.                                      final ODEStateAndDerivative globalPreviousState,
  86.                                      final ODEStateAndDerivative globalCurrentState,
  87.                                      final ODEStateAndDerivative softPreviousState,
  88.                                      final ODEStateAndDerivative softCurrentState,
  89.                                      final EquationsMapper mapper) {
  90.         super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
  91.     }

  93.     /** {@inheritDoc} */
  94.     @Override
  95.     protected DormandPrince54StateInterpolator create(final boolean newForward, final double[][] newYDotK,
  96.                                                       final ODEStateAndDerivative newGlobalPreviousState,
  97.                                                       final ODEStateAndDerivative newGlobalCurrentState,
  98.                                                       final ODEStateAndDerivative newSoftPreviousState,
  99.                                                       final ODEStateAndDerivative newSoftCurrentState,
  100.                                                       final EquationsMapper newMapper) {
  101.         return new DormandPrince54StateInterpolator(newForward, newYDotK,
  102.                                                     newGlobalPreviousState, newGlobalCurrentState,
  103.                                                     newSoftPreviousState, newSoftCurrentState,
  104.                                                     newMapper);
  105.     }

  107.     /** {@inheritDoc} */
  108.     @Override
  109.     protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
  110.                                                                            final double time, final double theta,
  111.                                                                            final double thetaH, final double oneMinusThetaH) {

  113.         // interpolate
  114.         final double eta = 1 - theta;
  115.         final double twoTheta = 2 * theta;
  116.         final double dot2 = 1 - twoTheta;
  117.         final double dot3 = theta * (2 - 3 * theta);
  118.         final double dot4 = twoTheta * (1 + theta * (twoTheta - 3));

  120.         final double[] interpolatedState;
  121.         final double[] interpolatedDerivatives;
  122.         if (getGlobalPreviousState() != null && theta <= 0.5) {
  123.             final double f1        = thetaH;
  124.             final double f2        = f1 * eta;
  125.             final double f3        = f2 * theta;
  126.             final double f4        = f3 * eta;
  127.             final double coeff0    = f1 * A70 - f2   * (A70 - 1) + f3   * (2 * A70 - 1) + f4   * D0;
  128.             final double coeff1    = 0;
  129.             final double coeff2    = f1 * A72 - f2   * A72       + f3   * (2 * A72)     + f4   * D2;
  130.             final double coeff3    = f1 * A73 - f2   * A73       + f3   * (2 * A73)     + f4   * D3;
  131.             final double coeff4    = f1 * A74 - f2   * A74       + f3   * (2 * A74)     + f4   * D4;
  132.             final double coeff5    = f1 * A75 - f2   * A75       + f3   * (2 * A75)     + f4   * D5;
  133.             final double coeff6    = f4 * D6 - f3;
  134.             final double coeffDot0 =      A70 - dot2 * (A70 - 1) + dot3 * (2 * A70 - 1) + dot4 * D0;
  135.             final double coeffDot1 = 0;
  136.             final double coeffDot2 =      A72 - dot2 * A72       + dot3 * (2 * A72)     + dot4 * D2;
  137.             final double coeffDot3 =      A73 - dot2 * A73       + dot3 * (2 * A73)     + dot4 * D3;
  138.             final double coeffDot4 =      A74 - dot2 * A74       + dot3 * (2 * A74)     + dot4 * D4;
  139.             final double coeffDot5 =      A75 - dot2 * A75       + dot3 * (2 * A75)     + dot4 * D5;
  140.             final double coeffDot6 = dot4 * D6 - dot3;
  141.             interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  142.                                                                      coeff4, coeff5, coeff6);
  143.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  144.                                                                   coeffDot4, coeffDot5, coeffDot6);
  145.         } else {
  146.             final double f1        = -oneMinusThetaH;
  147.             final double f2        = oneMinusThetaH * theta;
  148.             final double f3        = f2 * theta;
  149.             final double f4        = f3 * eta;
  150.             final double coeff0    = f1 * A70 - f2   * (A70 - 1) + f3   * (2 * A70 - 1) + f4   * D0;
  151.             final double coeff1    = 0;
  152.             final double coeff2    = f1 * A72 - f2   * A72       + f3   * (2 * A72)     + f4   * D2;
  153.             final double coeff3    = f1 * A73 - f2   * A73       + f3   * (2 * A73)     + f4   * D3;
  154.             final double coeff4    = f1 * A74 - f2   * A74       + f3   * (2 * A74)     + f4   * D4;
  155.             final double coeff5    = f1 * A75 - f2   * A75       + f3   * (2 * A75)     + f4   * D5;
  156.             final double coeff6    = f4 * D6 - f3;
  157.             final double coeffDot0 =      A70 - dot2 * (A70 - 1) + dot3 * (2 * A70 - 1) + dot4 * D0;
  158.             final double coeffDot1 = 0;
  159.             final double coeffDot2 =      A72 - dot2 * A72       + dot3 * (2 * A72)     + dot4 * D2;
  160.             final double coeffDot3 =      A73 - dot2 * A73       + dot3 * (2 * A73)     + dot4 * D3;
  161.             final double coeffDot4 =      A74 - dot2 * A74       + dot3 * (2 * A74)     + dot4 * D4;
  162.             final double coeffDot5 =      A75 - dot2 * A75       + dot3 * (2 * A75)     + dot4 * D5;
  163.             final double coeffDot6 = dot4 * D6 - dot3;
  164.             interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  165.                                                                     coeff4, coeff5, coeff6);
  166.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  167.                                                                   coeffDot4, coeffDot5, coeffDot6);
  168.         }

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

  172.     }

  174. }
Created with JaCoCo 0.8.12.202403310830