DormandPrince54FieldStepInterpolator.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.  */

  17. package org.apache.commons.math4.legacy.ode.nonstiff;

  18. import org.apache.commons.math4.legacy.core.Field;
  19. import org.apache.commons.math4.legacy.core.RealFieldElement;
  20. import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
  21. import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;

  22. /**
  23.  * This class represents an interpolator over the last step during an
  24.  * ODE integration for the 5(4) Dormand-Prince integrator.
  25.  *
  26.  * @see DormandPrince54Integrator
  27.  *
  28.  * @param <T> the type of the field elements
  29.  * @since 3.6
  30.  */

  31. class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
  32.       extends RungeKuttaFieldStepInterpolator<T> {

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

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

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

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

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

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

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

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

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

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

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

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

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

  57.     /** Simple constructor.
  58.      * @param field field to which the time and state vector elements belong
  59.      * @param forward integration direction indicator
  60.      * @param yDotK slopes at the intermediate points
  61.      * @param globalPreviousState start of the global step
  62.      * @param globalCurrentState end of the global step
  63.      * @param softPreviousState start of the restricted step
  64.      * @param softCurrentState end of the restricted step
  65.      * @param mapper equations mapper for the all equations
  66.      */
  67.     DormandPrince54FieldStepInterpolator(final Field<T> field, final boolean forward,
  68.                                          final T[][] yDotK,
  69.                                          final FieldODEStateAndDerivative<T> globalPreviousState,
  70.                                          final FieldODEStateAndDerivative<T> globalCurrentState,
  71.                                          final FieldODEStateAndDerivative<T> softPreviousState,
  72.                                          final FieldODEStateAndDerivative<T> softCurrentState,
  73.                                          final FieldEquationsMapper<T> mapper) {
  74.         super(field, forward, yDotK,
  75.               globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
  76.               mapper);
  77.         final T one = field.getOne();
  78.         a70 = one.multiply(   35.0).divide( 384.0);
  79.         a72 = one.multiply(  500.0).divide(1113.0);
  80.         a73 = one.multiply(  125.0).divide( 192.0);
  81.         a74 = one.multiply(-2187.0).divide(6784.0);
  82.         a75 = one.multiply(   11.0).divide(  84.0);
  83.         d0  = one.multiply(-12715105075.0).divide( 11282082432.0);
  84.         d2  = one.multiply( 87487479700.0).divide( 32700410799.0);
  85.         d3  = one.multiply(-10690763975.0).divide(  1880347072.0);
  86.         d4  = one.multiply(701980252875.0).divide(199316789632.0);
  87.         d5  = one.multiply( -1453857185.0).divide(   822651844.0);
  88.         d6  = one.multiply(    69997945.0).divide(    29380423.0);
  89.     }

  90.     /** {@inheritDoc} */
  91.     @Override
  92.     protected DormandPrince54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
  93.                                                                  final FieldODEStateAndDerivative<T> newGlobalPreviousState,
  94.                                                                  final FieldODEStateAndDerivative<T> newGlobalCurrentState,
  95.                                                                  final FieldODEStateAndDerivative<T> newSoftPreviousState,
  96.                                                                  final FieldODEStateAndDerivative<T> newSoftCurrentState,
  97.                                                                  final FieldEquationsMapper<T> newMapper) {
  98.         return new DormandPrince54FieldStepInterpolator<>(newField, newForward, newYDotK,
  99.                                                            newGlobalPreviousState, newGlobalCurrentState,
  100.                                                            newSoftPreviousState, newSoftCurrentState,
  101.                                                            newMapper);
  102.     }
  103.     /** {@inheritDoc} */
  104.     @SuppressWarnings("unchecked")
  105.     @Override
  106.     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
  107.                                                                                    final T time, final T theta,
  108.                                                                                    final T thetaH, final T oneMinusThetaH) {

  109.         // interpolate
  110.         final T one      = time.getField().getOne();
  111.         final T eta      = one.subtract(theta);
  112.         final T twoTheta = theta.multiply(2);
  113.         final T dot2     = one.subtract(twoTheta);
  114.         final T dot3     = theta.multiply(theta.multiply(-3).add(2));
  115.         final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
  116.         final T[] interpolatedState;
  117.         final T[] interpolatedDerivatives;
  118.         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
  119.             final T f1        = thetaH;
  120.             final T f2        = f1.multiply(eta);
  121.             final T f3        = f2.multiply(theta);
  122.             final T f4        = f3.multiply(eta);
  123.             final T coeff0    = f1.multiply(a70).
  124.                                 subtract(f2.multiply(a70.subtract(1))).
  125.                                 add(f3.multiply(a70.multiply(2).subtract(1))).
  126.                                 add(f4.multiply(d0));
  127.             final T coeff1    = time.getField().getZero();
  128.             final T coeff2    = f1.multiply(a72).
  129.                                 subtract(f2.multiply(a72)).
  130.                                 add(f3.multiply(a72.multiply(2))).
  131.                                 add(f4.multiply(d2));
  132.             final T coeff3    = f1.multiply(a73).
  133.                                 subtract(f2.multiply(a73)).
  134.                                 add(f3.multiply(a73.multiply(2))).
  135.                                 add(f4.multiply(d3));
  136.             final T coeff4    = f1.multiply(a74).
  137.                                 subtract(f2.multiply(a74)).
  138.                                 add(f3.multiply(a74.multiply(2))).
  139.                                 add(f4.multiply(d4));
  140.             final T coeff5    = f1.multiply(a75).
  141.                                 subtract(f2.multiply(a75)).
  142.                                 add(f3.multiply(a75.multiply(2))).
  143.                                 add(f4.multiply(d5));
  144.             final T coeff6    = f4.multiply(d6).subtract(f3);
  145.             final T coeffDot0 = a70.
  146.                                 subtract(dot2.multiply(a70.subtract(1))).
  147.                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
  148.                                 add(dot4.multiply(d0));
  149.             final T coeffDot1 = time.getField().getZero();
  150.             final T coeffDot2 = a72.
  151.                                 subtract(dot2.multiply(a72)).
  152.                                 add(dot3.multiply(a72.multiply(2))).
  153.                                 add(dot4.multiply(d2));
  154.             final T coeffDot3 = a73.
  155.                                 subtract(dot2.multiply(a73)).
  156.                                 add(dot3.multiply(a73.multiply(2))).
  157.                                 add(dot4.multiply(d3));
  158.             final T coeffDot4 = a74.
  159.                                 subtract(dot2.multiply(a74)).
  160.                                 add(dot3.multiply(a74.multiply(2))).
  161.                                 add(dot4.multiply(d4));
  162.             final T coeffDot5 = a75.
  163.                                 subtract(dot2.multiply(a75)).
  164.                                 add(dot3.multiply(a75.multiply(2))).
  165.                                 add(dot4.multiply(d5));
  166.             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
  167.             interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  168.                                                                      coeff4, coeff5, coeff6);
  169.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  170.                                                                   coeffDot4, coeffDot5, coeffDot6);
  171.         } else {
  172.             final T f1        = oneMinusThetaH.negate();
  173.             final T f2        = oneMinusThetaH.multiply(theta);
  174.             final T f3        = f2.multiply(theta);
  175.             final T f4        = f3.multiply(eta);
  176.             final T coeff0    = f1.multiply(a70).
  177.                                 subtract(f2.multiply(a70.subtract(1))).
  178.                                 add(f3.multiply(a70.multiply(2).subtract(1))).
  179.                                 add(f4.multiply(d0));
  180.             final T coeff1    = time.getField().getZero();
  181.             final T coeff2    = f1.multiply(a72).
  182.                                 subtract(f2.multiply(a72)).
  183.                                 add(f3.multiply(a72.multiply(2))).
  184.                                 add(f4.multiply(d2));
  185.             final T coeff3    = f1.multiply(a73).
  186.                                 subtract(f2.multiply(a73)).
  187.                                 add(f3.multiply(a73.multiply(2))).
  188.                                 add(f4.multiply(d3));
  189.             final T coeff4    = f1.multiply(a74).
  190.                                 subtract(f2.multiply(a74)).
  191.                                 add(f3.multiply(a74.multiply(2))).
  192.                                 add(f4.multiply(d4));
  193.             final T coeff5    = f1.multiply(a75).
  194.                                 subtract(f2.multiply(a75)).
  195.                                 add(f3.multiply(a75.multiply(2))).
  196.                                 add(f4.multiply(d5));
  197.             final T coeff6    = f4.multiply(d6).subtract(f3);
  198.             final T coeffDot0 = a70.
  199.                                 subtract(dot2.multiply(a70.subtract(1))).
  200.                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
  201.                                 add(dot4.multiply(d0));
  202.             final T coeffDot1 = time.getField().getZero();
  203.             final T coeffDot2 = a72.
  204.                                 subtract(dot2.multiply(a72)).
  205.                                 add(dot3.multiply(a72.multiply(2))).
  206.                                 add(dot4.multiply(d2));
  207.             final T coeffDot3 = a73.
  208.                                 subtract(dot2.multiply(a73)).
  209.                                 add(dot3.multiply(a73.multiply(2))).
  210.                                 add(dot4.multiply(d3));
  211.             final T coeffDot4 = a74.
  212.                                 subtract(dot2.multiply(a74)).
  213.                                 add(dot3.multiply(a74.multiply(2))).
  214.                                 add(dot4.multiply(d4));
  215.             final T coeffDot5 = a75.
  216.                                 subtract(dot2.multiply(a75)).
  217.                                 add(dot3.multiply(a75.multiply(2))).
  218.                                 add(dot4.multiply(d5));
  219.             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
  220.             interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
  221.                                                                     coeff4, coeff5, coeff6);
  222.             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
  223.                                                                   coeffDot4, coeffDot5, coeffDot6);
  224.         }
  225.         return new FieldODEStateAndDerivative<>(time, interpolatedState, interpolatedDerivatives);
  226.     }
  227. }