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

  17. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */

  21. package org.hipparchus.ode.nonstiff;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.ode.FieldEquationsMapper;
  25. import org.hipparchus.ode.FieldODEStateAndDerivative;

  26. /**
  27.  * This class represents an interpolator over the last step during an
  28.  * ODE integration for the 6th order Luther integrator.
  29.  *
  30.  * <p>This interpolator computes dense output inside the last
  31.  * step computed. The interpolation equation is consistent with the
  32.  * integration scheme.</p>
  33.  *
  34.  * @see LutherFieldIntegrator
  35.  * @param <T> the type of the field elements
  36.  */

  37. class LutherFieldStateInterpolator<T extends CalculusFieldElement<T>>
  38.     extends RungeKuttaFieldStateInterpolator<T> {

  39.     /** -49 - 49 q. */
  40.     private final T c5a;

  41.     /** 392 + 287 q. */
  42.     private final T c5b;

  43.     /** -637 - 357 q. */
  44.     private final T c5c;

  45.     /** 833 + 343 q. */
  46.     private final T c5d;

  47.     /** -49 + 49 q. */
  48.     private final T c6a;

  49.     /** -392 - 287 q. */
  50.     private final T c6b;

  51.     /** -637 + 357 q. */
  52.     private final T c6c;

  53.     /** 833 - 343 q. */
  54.     private final T c6d;

  55.     /** 49 + 49 q. */
  56.     private final T d5a;

  57.     /** -1372 - 847 q. */
  58.     private final T d5b;

  59.     /** 2254 + 1029 q */
  60.     private final T d5c;

  61.     /** 49 - 49 q. */
  62.     private final T d6a;

  63.     /** -1372 + 847 q. */
  64.     private final T d6b;

  65.     /** 2254 - 1029 q */
  66.     private final T d6c;

  67.     /** Simple constructor.
  68.      * @param field field to which the time and state vector elements belong
  69.      * @param forward integration direction indicator
  70.      * @param yDotK slopes at the intermediate points
  71.      * @param globalPreviousState start of the global step
  72.      * @param globalCurrentState end of the global step
  73.      * @param softPreviousState start of the restricted step
  74.      * @param softCurrentState end of the restricted step
  75.      * @param mapper equations mapper for the all equations
  76.      */
  77.     LutherFieldStateInterpolator(final Field<T> field, final boolean forward,
  78.                                  final T[][] yDotK,
  79.                                  final FieldODEStateAndDerivative<T> globalPreviousState,
  80.                                  final FieldODEStateAndDerivative<T> globalCurrentState,
  81.                                  final FieldODEStateAndDerivative<T> softPreviousState,
  82.                                  final FieldODEStateAndDerivative<T> softCurrentState,
  83.                                  final FieldEquationsMapper<T> mapper) {
  84.         super(field, forward, yDotK,
  85.               globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
  86.               mapper);
  87.         final T q = field.getZero().add(21).sqrt();
  88.         c5a = q.multiply(  -49).add(  -49);
  89.         c5b = q.multiply(  287).add(  392);
  90.         c5c = q.multiply( -357).add( -637);
  91.         c5d = q.multiply(  343).add(  833);
  92.         c6a = q.multiply(   49).add(  -49);
  93.         c6b = q.multiply( -287).add(  392);
  94.         c6c = q.multiply(  357).add( -637);
  95.         c6d = q.multiply( -343).add(  833);
  96.         d5a = q.multiply(   49).add(   49);
  97.         d5b = q.multiply( -847).add(-1372);
  98.         d5c = q.multiply( 1029).add( 2254);
  99.         d6a = q.multiply(  -49).add(   49);
  100.         d6b = q.multiply(  847).add(-1372);
  101.         d6c = q.multiply(-1029).add( 2254);
  102.     }

  103.     /** {@inheritDoc} */
  104.     @Override
  105.     protected LutherFieldStateInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
  106.                                                      final FieldODEStateAndDerivative<T> newGlobalPreviousState,
  107.                                                      final FieldODEStateAndDerivative<T> newGlobalCurrentState,
  108.                                                      final FieldODEStateAndDerivative<T> newSoftPreviousState,
  109.                                                      final FieldODEStateAndDerivative<T> newSoftCurrentState,
  110.                                                      final FieldEquationsMapper<T> newMapper) {
  111.         return new LutherFieldStateInterpolator<T>(newField, newForward, newYDotK,
  112.                                                    newGlobalPreviousState, newGlobalCurrentState,
  113.                                                    newSoftPreviousState, newSoftCurrentState,
  114.                                                    newMapper);
  115.     }

  116.     /** {@inheritDoc} */
  117.     @SuppressWarnings("unchecked")
  118.     @Override
  119.     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
  120.                                                                                    final T time, final T theta,
  121.                                                                                    final T thetaH, final T oneMinusThetaH) {

  122.         // the coefficients below have been computed by solving the
  123.         // order conditions from a theorem from Butcher (1963), using
  124.         // the method explained in Folkmar Bornemann paper "Runge-Kutta
  125.         // Methods, Trees, and Maple", Center of Mathematical Sciences, Munich
  126.         // University of Technology, February 9, 2001
  127.         //<http://wwwzenger.informatik.tu-muenchen.de/selcuk/sjam012101.html>

  128.         // the method is implemented in the rkcheck tool
  129.         // <https://www.spaceroots.org/software/rkcheck/index.html>.
  130.         // Running it for order 5 gives the following order conditions
  131.         // for an interpolator:
  132.         // order 1 conditions
  133.         // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
  134.         // order 2 conditions
  135.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
  136.         // order 3 conditions
  137.         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{2}}{6}
  138.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
  139.         // order 4 conditions
  140.         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{3}}{24}
  141.         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{3}}{12}
  142.         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{3}}{8}
  143.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
  144.         // order 5 conditions
  145.         // \sum_{i=4}^{i=s}\left(b_{i} \sum_{j=3}^{j=i-1}{\left(a_{i,j} \sum_{k=2}^{k=j-1}{\left(a_{j,k} \sum_{l=1}^{l=k-1}{\left(a_{k,l} c_{l} \right)} \right)} \right)}\right) = \frac{\theta^{4}}{120}
  146.         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k}^{2} \right)} \right)}\right) = \frac{\theta^{4}}{60}
  147.         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} c_{j}\sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{40}
  148.         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{3} \right)}\right) = \frac{\theta^{4}}{20}
  149.         // \sum_{i=3}^{i=s}\left(b_{i} c_{i}\sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{30}
  150.         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{4}}{15}
  151.         // \sum_{i=2}^{i=s}\left(b_{i} \left(\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)} \right)^{2}\right) = \frac{\theta^{4}}{20}
  152.         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}^{2}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{4}}{10}
  153.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}

  154.         // The a_{j,k} and c_{k} are given by the integrator Butcher arrays. What remains to solve
  155.         // are the b_i for the interpolator. They are found by solving the above equations.
  156.         // For a given interpolator, some equations are redundant, so in our case when we select
  157.         // all equations from order 1 to 4, we still don't have enough independent equations
  158.         // to solve from b_1 to b_7. We need to also select one equation from order 5. Here,
  159.         // we selected the last equation. It appears this choice implied at least the last 3 equations
  160.         // are fulfilled, but some of the former ones are not, so the resulting interpolator is order 5.
  161.         // At the end, we get the b_i as polynomials in theta.

  162.         final T coeffDot1 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(   21        ).add( -47          )).add(   36         )).add( -54     /   5.0)).add(1);
  163.         final T coeffDot2 =  time.getField().getZero();
  164.         final T coeffDot3 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(  112        ).add(-608    /  3.0)).add(  320   / 3.0 )).add(-208    /  15.0));
  165.         final T coeffDot4 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567  /  5.0).add( 972    /  5.0)).add( -486   / 5.0 )).add( 324    /  25.0));
  166.         final T coeffDot5 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150)));
  167.         final T coeffDot6 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150)));
  168.         final T coeffDot7 =  theta.multiply(theta.multiply(theta.multiply(                                             3.0 ).add(   -3         )).add(   3   /   5.0));
  169.         final T[] interpolatedState;
  170.         final T[] interpolatedDerivatives;

  171.         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {

  172.             final T s         = thetaH;
  173.             final T coeff1    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(  21    /  5.0).add( -47    /  4.0)).add(   12         )).add( -27    /   5.0)).add(1));
  174.             final T coeff2    = time.getField().getZero();
  175.             final T coeff3    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112    /  5.0).add(-152    /  3.0)).add(  320   / 9.0 )).add(-104    /  15.0)));
  176.             final T coeff4    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567    / 25.0).add( 243    /  5.0)).add( -162   / 5.0 )).add( 162    /  25.0)));
  177.             final T coeff5    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))));
  178.             final T coeff6    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))));
  179.             final T coeff7    = s.multiply(theta.multiply(theta.multiply(theta.multiply(                                      3    /  4.0 ).add(   -1         )).add(   3    /  10.0)));
  180.             interpolatedState       = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
  181.             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
  182.         } else {

  183.             final T s         = oneMinusThetaH;
  184.             final T coeff1    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21   /   5.0).add(   151  /  20.0)).add( -89   /   20.0)).add(  19 /  20.0)).add(- 1 / 20.0));
  185.             final T coeff2    = time.getField().getZero();
  186.             final T coeff3    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112   /   5.0).add(   424  /  15.0)).add( -328  /   45.0)).add( -16 /  45.0)).add(-16 /  45.0));
  187.             final T coeff4    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567   /  25.0).add(  -648  /  25.0)).add(  162  /   25.0))));
  188.             final T coeff5    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
  189.             final T coeff6    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
  190.             final T coeff7    = s.multiply(               theta.multiply(theta.multiply(theta.multiply(                        -3  /   4.0 ).add(   1   /    4.0)).add(  -1 /  20.0)).add( -1 /  20.0));
  191.             interpolatedState       = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
  192.             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
  193.         }

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

  195.     }

  196. }