DormandPrince54FieldStateInterpolator.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.ode.FieldEquationsMapper;
  28. import org.hipparchus.ode.FieldODEStateAndDerivative;
  29. import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;

  31. /**
  32.  * This class represents an interpolator over the last step during an
  33.  * ODE integration for the 5(4) Dormand-Prince integrator.
  34.  *
  35.  * @see DormandPrince54Integrator
  36.  *
  37.  * @param <T> the type of the field elements
  38.  */

  40. public class DormandPrince54FieldStateInterpolator<T extends CalculusFieldElement<T>>
  41.     extends RungeKuttaFieldStateInterpolator<T> {

  43.     /** Last row of the Butcher-array internal weights, element 0. */
  44.     private final T a70;

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

  48.     /** Last row of the Butcher-array internal weights, element 2. */
  49.     private final T a72;

  51.     /** Last row of the Butcher-array internal weights, element 3. */
  52.     private final T a73;

  54.     /** Last row of the Butcher-array internal weights, element 4. */
  55.     private final T a74;

  57.     /** Last row of the Butcher-array internal weights, element 5. */
  58.     private final T a75;

  60.     /** Shampine (1986) Dense output, element 0. */
  61.     private final T d0;

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

  65.     /** Shampine (1986) Dense output, element 2. */
  66.     private final T d2;

  68.     /** Shampine (1986) Dense output, element 3. */
  69.     private final T d3;

  71.     /** Shampine (1986) Dense output, element 4. */
  72.     private final T d4;

  74.     /** Shampine (1986) Dense output, element 5. */
  75.     private final T d5;

  77.     /** Shampine (1986) Dense output, element 6. */
  78.     private final T d6;

  80.     /** Simple constructor.
  81.      * @param field field to which the time and state vector elements belong
  82.      * @param forward integration direction indicator
  83.      * @param yDotK slopes at the intermediate points
  84.      * @param globalPreviousState start of the global step
  85.      * @param globalCurrentState end of the global step
  86.      * @param softPreviousState start of the restricted step
  87.      * @param softCurrentState end of the restricted step
  88.      * @param mapper equations mapper for the all equations
  89.      */
  90.     public DormandPrince54FieldStateInterpolator(final Field<T> field, final boolean forward,
  91.                                                  final T[][] yDotK,
  92.                                                  final FieldODEStateAndDerivative<T> globalPreviousState,
  93.                                                  final FieldODEStateAndDerivative<T> globalCurrentState,
  94.                                                  final FieldODEStateAndDerivative<T> softPreviousState,
  95.                                                  final FieldODEStateAndDerivative<T> softCurrentState,
  96.                                                  final FieldEquationsMapper<T> mapper) {
  97.         super(field, forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
  98.                 mapper);
  99.         final T one = field.getOne();
  100.         a70 = one.newInstance(   35.0 /  384.0);
  101.         a72 = one.newInstance(  500.0 / 1113.0);
  102.         a73 = one.newInstance(  125.0 /  192.0);
  103.         a74 = one.newInstance(-2187.0 / 6784.0);
  104.         a75 = one.newInstance(   11.0 /   84.0);
  105.         d0  = one.newInstance(-12715105075.0 /  11282082432.0);
  106.         d2  = one.newInstance( 87487479700.0 /  32700410799.0);
  107.         d3  = one.newInstance(-10690763975.0 /   1880347072.0);
  108.         d4  = one.newInstance(701980252875.0 / 199316789632.0);
  109.         d5  = one.newInstance( -1453857185.0 /    822651844.0);
  110.         d6  = one.newInstance(    69997945.0 /     29380423.0);
  111.     }

  113.     /** {@inheritDoc} */
  114.     @Override
  115.     protected DormandPrince54FieldStateInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
  116.                                                               final FieldODEStateAndDerivative<T> newGlobalPreviousState,
  117.                                                               final FieldODEStateAndDerivative<T> newGlobalCurrentState,
  118.                                                               final FieldODEStateAndDerivative<T> newSoftPreviousState,
  119.                                                               final FieldODEStateAndDerivative<T> newSoftCurrentState,
  120.                                                               final FieldEquationsMapper<T> newMapper) {
  121.         return new DormandPrince54FieldStateInterpolator<>(newField, newForward, newYDotK,
  122.                                                             newGlobalPreviousState, newGlobalCurrentState,
  123.                                                             newSoftPreviousState, newSoftCurrentState,
  124.                                                             newMapper);
  125.     }
  126.     /** {@inheritDoc} */
  127.     @SuppressWarnings("unchecked")
  128.     @Override
  129.     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
  130.                                                                                    final T time, final T theta,
  131.                                                                                    final T thetaH, final T oneMinusThetaH) {

  133.         // interpolate
  134.         final T one      = time.getField().getOne();
  135.         final T eta      = one.subtract(theta);
  136.         final T twoTheta = theta.multiply(2);
  137.         final T dot2     = one.subtract(twoTheta);
  138.         final T dot3     = theta.multiply(theta.multiply(-3).add(2));
  139.         final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
  140.         final T[] interpolatedState;
  141.         final T[] interpolatedDerivatives;
  142.         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
  143.             final T f1        = thetaH;
  144.             final T f2        = f1.multiply(eta);
  145.             final T f3        = f2.multiply(theta);
  146.             final T f4        = f3.multiply(eta);
  147.             final T coeff0    = f1.multiply(a70).
  148.                                 subtract(f2.multiply(a70.subtract(1))).
  149.                                 add(f3.multiply(a70.multiply(2).subtract(1))).
  150.                                 add(f4.multiply(d0));
  151.             final T coeff1    = time.getField().getZero();
  152.             final T coeff2    = f1.multiply(a72).
  153.                                 subtract(f2.multiply(a72)).
  154.                                 add(f3.multiply(a72.multiply(2))).
  155.                                 add(f4.multiply(d2));
  156.             final T coeff3    = f1.multiply(a73).
  157.                                 subtract(f2.multiply(a73)).
  158.                                 add(f3.multiply(a73.multiply(2))).
  159.                                 add(f4.multiply(d3));
  160.             final T coeff4    = f1.multiply(a74).
  161.                                 subtract(f2.multiply(a74)).
  162.                                 add(f3.multiply(a74.multiply(2))).
  163.                                 add(f4.multiply(d4));
  164.             final T coeff5    = f1.multiply(a75).
  165.                                 subtract(f2.multiply(a75)).
  166.                                 add(f3.multiply(a75.multiply(2))).
  167.                                 add(f4.multiply(d5));
  168.             final T coeff6    = f4.multiply(d6).subtract(f3);
  169.             final T coeffDot0 = a70.
  170.                                 subtract(dot2.multiply(a70.subtract(1))).
  171.                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
  172.                                 add(dot4.multiply(d0));
  173.             final T coeffDot1 = time.getField().getZero();
  174.             final T coeffDot2 = a72.
  175.                                 subtract(dot2.multiply(a72)).
  176.                                 add(dot3.multiply(a72.multiply(2))).
  177.                                 add(dot4.multiply(d2));
  178.             final T coeffDot3 = a73.
  179.                                 subtract(dot2.multiply(a73)).
  180.                                 add(dot3.multiply(a73.multiply(2))).
  181.                                 add(dot4.multiply(d3));
  182.             final T coeffDot4 = a74.
  183.                                 subtract(dot2.multiply(a74)).
  184.                                 add(dot3.multiply(a74.multiply(2))).
  185.                                 add(dot4.multiply(d4));
  186.             final T coeffDot5 = a75.
  187.                                 subtract(dot2.multiply(a75)).
  188.                                 add(dot3.multiply(a75.multiply(2))).
  189.                                 add(dot4.multiply(d5));
  190.             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
  191.             interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  192.                                                                      coeff4, coeff5, coeff6);
  193.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  194.                                                                   coeffDot4, coeffDot5, coeffDot6);
  195.         } else {
  196.             final T f1        = oneMinusThetaH.negate();
  197.             final T f2        = oneMinusThetaH.multiply(theta);
  198.             final T f3        = f2.multiply(theta);
  199.             final T f4        = f3.multiply(eta);
  200.             final T coeff0    = f1.multiply(a70).
  201.                                 subtract(f2.multiply(a70.subtract(1))).
  202.                                 add(f3.multiply(a70.multiply(2).subtract(1))).
  203.                                 add(f4.multiply(d0));
  204.             final T coeff1    = time.getField().getZero();
  205.             final T coeff2    = f1.multiply(a72).
  206.                                 subtract(f2.multiply(a72)).
  207.                                 add(f3.multiply(a72.multiply(2))).
  208.                                 add(f4.multiply(d2));
  209.             final T coeff3    = f1.multiply(a73).
  210.                                 subtract(f2.multiply(a73)).
  211.                                 add(f3.multiply(a73.multiply(2))).
  212.                                 add(f4.multiply(d3));
  213.             final T coeff4    = f1.multiply(a74).
  214.                                 subtract(f2.multiply(a74)).
  215.                                 add(f3.multiply(a74.multiply(2))).
  216.                                 add(f4.multiply(d4));
  217.             final T coeff5    = f1.multiply(a75).
  218.                                 subtract(f2.multiply(a75)).
  219.                                 add(f3.multiply(a75.multiply(2))).
  220.                                 add(f4.multiply(d5));
  221.             final T coeff6    = f4.multiply(d6).subtract(f3);
  222.             final T coeffDot0 = a70.
  223.                                 subtract(dot2.multiply(a70.subtract(1))).
  224.                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
  225.                                 add(dot4.multiply(d0));
  226.             final T coeffDot1 = time.getField().getZero();
  227.             final T coeffDot2 = a72.
  228.                                 subtract(dot2.multiply(a72)).
  229.                                 add(dot3.multiply(a72.multiply(2))).
  230.                                 add(dot4.multiply(d2));
  231.             final T coeffDot3 = a73.
  232.                                 subtract(dot2.multiply(a73)).
  233.                                 add(dot3.multiply(a73.multiply(2))).
  234.                                 add(dot4.multiply(d3));
  235.             final T coeffDot4 = a74.
  236.                                 subtract(dot2.multiply(a74)).
  237.                                 add(dot3.multiply(a74.multiply(2))).
  238.                                 add(dot4.multiply(d4));
  239.             final T coeffDot5 = a75.
  240.                                 subtract(dot2.multiply(a75)).
  241.                                 add(dot3.multiply(a75.multiply(2))).
  242.                                 add(dot4.multiply(d5));
  243.             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
  244.             interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  245.                                                                     coeff4, coeff5, coeff6);
  246.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  247.                                                                   coeffDot4, coeffDot5, coeffDot6);
  248.         }
  249.         return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);

  251.     }

  253. }
Created with JaCoCo 0.8.12.202403310830