EmbeddedRungeKuttaIntegrator.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;

  20. import org.hipparchus.exception.MathIllegalArgumentException;
  21. import org.hipparchus.exception.MathIllegalStateException;
  22. import org.hipparchus.ode.EquationsMapper;
  23. import org.hipparchus.ode.ExpandableODE;
  24. import org.hipparchus.ode.LocalizedODEFormats;
  25. import org.hipparchus.ode.ODEState;
  26. import org.hipparchus.ode.ODEStateAndDerivative;
  27. import org.hipparchus.ode.nonstiff.interpolators.RungeKuttaStateInterpolator;
  28. import org.hipparchus.util.FastMath;

  30. /**
  31.  * This class implements the common part of all embedded Runge-Kutta
  32.  * integrators for Ordinary Differential Equations.
  33.  *
  34.  * <p>These methods are embedded explicit Runge-Kutta methods with two
  35.  * sets of coefficients allowing to estimate the error, their Butcher
  36.  * arrays are as follows :</p>
  37.  * <pre>
  38.  *    0  |
  39.  *   c2  | a21
  40.  *   c3  | a31  a32
  41.  *   ... |        ...
  42.  *   cs  | as1  as2  ...  ass-1
  43.  *       |--------------------------
  44.  *       |  b1   b2  ...   bs-1  bs
  45.  *       |  b'1  b'2 ...   b's-1 b's
  46.  * </pre>
  47.  *
  48.  * <p>In fact, we rather use the array defined by ej = bj - b'j to
  49.  * compute directly the error rather than computing two estimates and
  50.  * then comparing them.</p>
  51.  *
  52.  * <p>Some methods are qualified as <i>fsal</i> (first same as last)
  53.  * methods. This means the last evaluation of the derivatives in one
  54.  * step is the same as the first in the next step. Then, this
  55.  * evaluation can be reused from one step to the next one and the cost
  56.  * of such a method is really s-1 evaluations despite the method still
  57.  * has s stages. This behaviour is true only for successful steps, if
  58.  * the step is rejected after the error estimation phase, no
  59.  * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
  60.  * asi = bi for all i.</p>
  61.  *
  62.  */

  64. public abstract class EmbeddedRungeKuttaIntegrator
  65.     extends AdaptiveStepsizeIntegrator
  66.     implements ExplicitRungeKuttaIntegrator {

  68.     /** Index of the pre-computed derivative for <i>fsal</i> methods. */
  69.     private final int fsal;

  71.     /** Time steps from Butcher array (without the first zero). */
  72.     private final double[] c;

  74.     /** Internal weights from Butcher array (without the first empty row). */
  75.     private final double[][] a;

  77.     /** External weights for the high order method from Butcher array. */
  78.     private final double[] b;

  80.     /** Stepsize control exponent. */
  81.     private final double exp;

  83.     /** Safety factor for stepsize control. */
  84.     private double safety;

  86.     /** Minimal reduction factor for stepsize control. */
  87.     private double minReduction;

  89.     /** Maximal growth factor for stepsize control. */
  90.     private double maxGrowth;

  92.     /** Build a Runge-Kutta integrator with the given Butcher array.
  93.      * @param name name of the method
  94.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  95.      * or -1 if method is not <i>fsal</i>
  96.      * @param minStep minimal step (sign is irrelevant, regardless of
  97.      * integration direction, forward or backward), the last step can
  98.      * be smaller than this
  99.      * @param maxStep maximal step (sign is irrelevant, regardless of
  100.      * integration direction, forward or backward), the last step can
  101.      * be smaller than this
  102.      * @param scalAbsoluteTolerance allowed absolute error
  103.      * @param scalRelativeTolerance allowed relative error
  104.      */
  105.     protected EmbeddedRungeKuttaIntegrator(final String name, final int fsal,
  106.                                            final double minStep, final double maxStep,
  107.                                            final double scalAbsoluteTolerance,
  108.                                            final double scalRelativeTolerance) {

  110.         super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

  112.         this.fsal = fsal;
  113.         this.c    = getC();
  114.         this.a    = getA();
  115.         this.b    = getB();

  117.         exp = -1.0 / getOrder();

  119.         // set the default values of the algorithm control parameters
  120.         setSafety(0.9);
  121.         setMinReduction(0.2);
  122.         setMaxGrowth(10.0);

  124.     }

  126.     /** Build a Runge-Kutta integrator with the given Butcher array.
  127.      * @param name name of the method
  128.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  129.      * or -1 if method is not <i>fsal</i>
  130.      * @param minStep minimal step (must be positive even for backward
  131.      * integration), the last step can be smaller than this
  132.      * @param maxStep maximal step (must be positive even for backward
  133.      * integration)
  134.      * @param vecAbsoluteTolerance allowed absolute error
  135.      * @param vecRelativeTolerance allowed relative error
  136.      */
  137.     protected EmbeddedRungeKuttaIntegrator(final String name, final int fsal,
  138.                                            final double   minStep, final double maxStep,
  139.                                            final double[] vecAbsoluteTolerance,
  140.                                            final double[] vecRelativeTolerance) {

  142.         super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);

  144.         this.fsal = fsal;
  145.         this.c    = getC();
  146.         this.a    = getA();
  147.         this.b    = getB();

  149.         exp = -1.0 / getOrder();

  151.         // set the default values of the algorithm control parameters
  152.         setSafety(0.9);
  153.         setMinReduction(0.2);
  154.         setMaxGrowth(10.0);

  156.     }

  158.     /** Create an interpolator.
  159.      * @param forward integration direction indicator
  160.      * @param yDotK slopes at the intermediate points
  161.      * @param globalPreviousState start of the global step
  162.      * @param globalCurrentState end of the global step
  163.      * @param mapper equations mapper for the all equations
  164.      * @return external weights for the high order method from Butcher array
  165.      */
  166.     protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
  167.                                                                       ODEStateAndDerivative globalPreviousState,
  168.                                                                       ODEStateAndDerivative globalCurrentState,
  169.                                                                       EquationsMapper mapper);
  170.     /** Get the order of the method.
  171.      * @return order of the method
  172.      */
  173.     public abstract int getOrder();

  175.     /** Get the safety factor for stepsize control.
  176.      * @return safety factor
  177.      */
  178.     public double getSafety() {
  179.         return safety;
  180.     }

  182.     /** Set the safety factor for stepsize control.
  183.      * @param safety safety factor
  184.      */
  185.     public void setSafety(final double safety) {
  186.         this.safety = safety;
  187.     }

  189.     /** {@inheritDoc} */
  190.     @Override
  191.     public ODEStateAndDerivative integrate(final ExpandableODE equations,
  192.                                            final ODEState initialState, final double finalTime)
  193.         throws MathIllegalArgumentException, MathIllegalStateException {

  195.         sanityChecks(initialState, finalTime);
  196.         setStepStart(initIntegration(equations, initialState, finalTime));
  197.         final boolean forward = finalTime > initialState.getTime();

  199.         // create some internal working arrays
  200.         final int        stages  = c.length + 1;
  201.         final double[][] yDotK   = new double[stages][];
  202.         double[]   yTmp    = new double[equations.getMapper().getTotalDimension()];

  204.         // set up integration control objects
  205.         double  hNew      = 0;
  206.         boolean firstTime = true;

  208.         // main integration loop
  209.         setIsLastStep(false);
  210.         do {

  212.             // iterate over step size, ensuring local normalized error is smaller than 1
  213.             double error = 10;
  214.             while (error >= 1.0) {

  216.                 // first stage
  217.                 final double[] y = getStepStart().getCompleteState();
  218.                 yDotK[0] = getStepStart().getCompleteDerivative();

  220.                 if (firstTime) {
  221.                     final StepsizeHelper helper = getStepSizeHelper();
  222.                     final double[] scale = new double[helper.getMainSetDimension()];
  223.                     for (int i = 0; i < scale.length; ++i) {
  224.                         scale[i] = helper.getTolerance(i, FastMath.abs(y[i]));
  225.                     }
  226.                     hNew = initializeStep(forward, getOrder(), scale, getStepStart());
  227.                     firstTime = false;
  228.                 }

  230.                 setStepSize(hNew);
  231.                 if (forward) {
  232.                     if (getStepStart().getTime() + getStepSize() >= finalTime) {
  233.                         setStepSize(finalTime - getStepStart().getTime());
  234.                     }
  235.                 } else {
  236.                     if (getStepStart().getTime() + getStepSize() <= finalTime) {
  237.                         setStepSize(finalTime - getStepStart().getTime());
  238.                     }
  239.                 }

  241.                 // next stages
  242.                 ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(), y,
  243.                         getStepSize(), a, c, yDotK);
  244.                 yTmp = ExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);

  246.                 incrementEvaluations(stages - 1);

  248.                 // estimate the error at the end of the step
  249.                 error = estimateError(yDotK, y, yTmp, getStepSize());
  250.                 if (Double.isNaN(error)) {
  251.                     throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
  252.                                                         getStepStart().getTime() + getStepSize());
  253.                 }
  254.                 if (error >= 1.0) {
  255.                     // reject the step and attempt to reduce error by stepsize control
  256.                     final double factor =
  257.                                     FastMath.min(maxGrowth,
  258.                                                  FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
  259.                     hNew = getStepSizeHelper().filterStep(getStepSize() * factor, forward, false);
  260.                 }

  262.             }
  263.             final double   stepEnd = getStepStart().getTime() + getStepSize();
  264.             final double[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
  265.             final ODEStateAndDerivative stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);

  267.             // local error is small enough: accept the step, trigger events and step handlers
  268.             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()), finalTime));

  270.             if (!isLastStep()) {

  272.                 // stepsize control for next step
  273.                 final double factor =
  274.                                 FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
  275.                 final double  scaledH    = getStepSize() * factor;
  276.                 final double  nextT      = getStepStart().getTime() + scaledH;
  277.                 final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
  278.                 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);

  280.                 final double  filteredNextT      = getStepStart().getTime() + hNew;
  281.                 final boolean filteredNextIsLast = forward ? (filteredNextT >= finalTime) : (filteredNextT <= finalTime);
  282.                 if (filteredNextIsLast) {
  283.                     hNew = finalTime - getStepStart().getTime();
  284.                 }

  286.             }

  288.         } while (!isLastStep());

  290.         final ODEStateAndDerivative finalState = getStepStart();
  291.         resetInternalState();
  292.         return finalState;

  294.     }

  296.     /** Get the minimal reduction factor for stepsize control.
  297.      * @return minimal reduction factor
  298.      */
  299.     public double getMinReduction() {
  300.         return minReduction;
  301.     }

  303.     /** Set the minimal reduction factor for stepsize control.
  304.      * @param minReduction minimal reduction factor
  305.      */
  306.     public void setMinReduction(final double minReduction) {
  307.         this.minReduction = minReduction;
  308.     }

  310.     /** Get the maximal growth factor for stepsize control.
  311.      * @return maximal growth factor
  312.      */
  313.     public double getMaxGrowth() {
  314.         return maxGrowth;
  315.     }

  317.     /** Set the maximal growth factor for stepsize control.
  318.      * @param maxGrowth maximal growth factor
  319.      */
  320.     public void setMaxGrowth(final double maxGrowth) {
  321.         this.maxGrowth = maxGrowth;
  322.     }

  324.     /** Compute the error ratio.
  325.      * @param yDotK derivatives computed during the first stages
  326.      * @param y0 estimate of the step at the start of the step
  327.      * @param y1 estimate of the step at the end of the step
  328.      * @param h  current step
  329.      * @return error ratio, greater than 1 if step should be rejected
  330.      */
  331.     protected abstract double estimateError(double[][] yDotK,
  332.                                             double[] y0, double[] y1,
  333.                                             double h);

  335. }
Created with JaCoCo 0.8.12.202403310830