LutherStateInterpolator.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.LutherIntegrator;
  23. import org.hipparchus.util.FastMath;

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

  36. public class LutherStateInterpolator extends RungeKuttaStateInterpolator {

  38.     /** Serializable version identifier */
  39.     private static final long serialVersionUID = 20160328;

  41.     /** Square root. */
  42.     private static final double Q = FastMath.sqrt(21);

  44.     /** Simple constructor.
  45.      * @param forward integration direction indicator
  46.      * @param yDotK slopes at the intermediate points
  47.      * @param globalPreviousState start of the global step
  48.      * @param globalCurrentState end of the global step
  49.      * @param softPreviousState start of the restricted step
  50.      * @param softCurrentState end of the restricted step
  51.      * @param mapper equations mapper for the all equations
  52.      */
  53.     public LutherStateInterpolator(final boolean forward,
  54.                                    final double[][] yDotK,
  55.                                    final ODEStateAndDerivative globalPreviousState,
  56.                                    final ODEStateAndDerivative globalCurrentState,
  57.                                    final ODEStateAndDerivative softPreviousState,
  58.                                    final ODEStateAndDerivative softCurrentState,
  59.                                    final EquationsMapper mapper) {
  60.         super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
  61.     }

  63.     /** {@inheritDoc} */
  64.     @Override
  65.     protected LutherStateInterpolator create(final boolean newForward, final double[][] newYDotK,
  66.                                              final ODEStateAndDerivative newGlobalPreviousState,
  67.                                              final ODEStateAndDerivative newGlobalCurrentState,
  68.                                              final ODEStateAndDerivative newSoftPreviousState,
  69.                                              final ODEStateAndDerivative newSoftCurrentState,
  70.                                              final EquationsMapper newMapper) {
  71.         return new LutherStateInterpolator(newForward, newYDotK,
  72.                                            newGlobalPreviousState, newGlobalCurrentState,
  73.                                            newSoftPreviousState, newSoftCurrentState,
  74.                                            newMapper);
  75.     }

  77.     /** {@inheritDoc} */
  78.     @Override
  79.     protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
  80.                                                                            final double time, final double theta,
  81.                                                                            final double thetaH, final double oneMinusThetaH) {

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

  90.         // the method is implemented in the rkcheck tool
  91.         // <https://www.spaceroots.org/software/rkcheck/index.html>.
  92.         // Running it for order 5 gives the following order conditions
  93.         // for an interpolator:
  94.         // order 1 conditions
  95.         // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
  96.         // order 2 conditions
  97.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
  98.         // order 3 conditions
  99.         // \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}
  100.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
  101.         // order 4 conditions
  102.         // \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}
  103.         // \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}
  104.         // \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}
  105.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
  106.         // order 5 conditions
  107.         // \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}
  108.         // \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}
  109.         // \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}
  110.         // \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}
  111.         // \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}
  112.         // \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}
  113.         // \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}
  114.         // \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}
  115.         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}

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

  126.         final double[] interpolatedState;
  127.         final double[] interpolatedDerivatives;

  129.         final double coeffDot1 =  1 + theta * ( -54            /   5.0 + theta * (   36                   + theta * ( -47                   + theta *   21)));
  130.         final double coeffDot2 =  0;
  131.         final double coeffDot3 =      theta * (-208            /  15.0 + theta * (  320            / 3.0  + theta * (-608            /  3.0 + theta *  112)));
  132.         final double coeffDot4 =      theta * ( 324            /  25.0 + theta * ( -486            / 5.0  + theta * ( 972            /  5.0 + theta * -567           /  5.0)));
  133.         final double coeffDot5 =      theta * ((833 + 343 * Q) / 150.0 + theta * ((-637 - 357 * Q) / 30.0 + theta * ((392 + 287 * Q) / 15.0 + theta * (-49 - 49 * Q) /  5.0)));
  134.         final double coeffDot6 =      theta * ((833 - 343 * Q) / 150.0 + theta * ((-637 + 357 * Q) / 30.0 + theta * ((392 - 287 * Q) / 15.0 + theta * (-49 + 49 * Q) /  5.0)));
  135.         final double coeffDot7 =      theta * (   3            /   5.0 + theta * (   -3                   + theta *     3));

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

  139.             final double coeff1    =  1 + theta * ( -27            /   5.0 + theta * (   12                   + theta * ( -47            /  4.0 + theta *   21           /  5.0)));
  140.             final double coeff2    =  0;
  141.             final double coeff3    =      theta * (-104            /  15.0 + theta * (  320            / 9.0  + theta * (-152            /  3.0 + theta *  112           /  5.0)));
  142.             final double coeff4    =      theta * ( 162            /  25.0 + theta * ( -162            / 5.0  + theta * ( 243            /  5.0 + theta * -567           / 25.0)));
  143.             final double coeff5    =      theta * ((833 + 343 * Q) / 300.0 + theta * ((-637 - 357 * Q) / 90.0 + theta * ((392 + 287 * Q) / 60.0 + theta * (-49 - 49 * Q) / 25.0)));
  144.             final double coeff6    =      theta * ((833 - 343 * Q) / 300.0 + theta * ((-637 + 357 * Q) / 90.0 + theta * ((392 - 287 * Q) / 60.0 + theta * (-49 + 49 * Q) / 25.0)));
  145.             final double coeff7    =      theta * (   3            /  10.0 + theta * (   -1                   + theta * (   3            /  4.0)));
  146.             interpolatedState       = previousStateLinearCombination(thetaH * coeff1, thetaH * coeff2,
  147.                                                                      thetaH * coeff3, thetaH * coeff4,
  148.                                                                      thetaH * coeff5, thetaH * coeff6,
  149.                                                                      thetaH * coeff7);
  150.             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
  151.         } else {

  153.             final double coeff1    =  -1 /  20.0 + theta * (  19            /  20.0 + theta * (  -89             /  20.0  + theta * (   151            /  20.0 + theta *  -21           /   5.0)));
  154.             final double coeff2    =  0;
  155.             final double coeff3    = -16 /  45.0 + theta * ( -16            /  45.0 + theta * ( -328             /  45.0  + theta * (   424            /  15.0 + theta * -112           /   5.0)));
  156.             final double coeff4    =               theta * (                          theta * (  162             /  25.0  + theta * (  -648            /  25.0 + theta *  567           /  25.0)));
  157.             final double coeff5    = -49 / 180.0 + theta * ( -49            / 180.0 + theta * ((2254 + 1029 * Q) / 900.0  + theta * ((-1372 - 847 * Q) / 300.0 + theta * ( 49 + 49 * Q) /  25.0)));
  158.             final double coeff6    = -49 / 180.0 + theta * ( -49            / 180.0 + theta * ((2254 - 1029 * Q) / 900.0  + theta * ((-1372 + 847 * Q) / 300.0 + theta * ( 49 - 49 * Q) /  25.0)));
  159.             final double coeff7    =  -1 /  20.0 + theta * (  -1            /  20.0 + theta * (    1             /   4.0  + theta * (    -3            /   4.0)));
  160.             interpolatedState       = currentStateLinearCombination(oneMinusThetaH * coeff1, oneMinusThetaH * coeff2,
  161.                                                                     oneMinusThetaH * coeff3, oneMinusThetaH * coeff4,
  162.                                                                     oneMinusThetaH * coeff5, oneMinusThetaH * coeff6,
  163.                                                                     oneMinusThetaH * coeff7);
  164.             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
  165.         }

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

  169.     }

  171. }
Created with JaCoCo 0.8.12.202403310830