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

  25. import org.hipparchus.CalculusFieldElement;
  26. import org.hipparchus.Field;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.exception.MathIllegalStateException;
  29. import org.hipparchus.ode.FieldEquationsMapper;
  30. import org.hipparchus.ode.FieldExpandableODE;
  31. import org.hipparchus.ode.FieldODEState;
  32. import org.hipparchus.ode.FieldODEStateAndDerivative;
  33. import org.hipparchus.ode.nonstiff.interpolators.RungeKuttaFieldStateInterpolator;
  34. import org.hipparchus.util.FastMath;
  35. import org.hipparchus.util.MathArrays;
  36. import org.hipparchus.util.MathUtils;

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

  73. public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
  74.     extends AdaptiveStepsizeFieldIntegrator<T>
  75.     implements FieldExplicitRungeKuttaIntegrator<T> {

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

  80.     /** Time steps from Butcher array (without the first zero). */
  81.     private final T[] c;

  83.     /** Internal weights from Butcher array (without the first empty row). */
  84.     private final T[][] a;

  86.     /** External weights for the high order method from Butcher array. */
  87.     private final T[] b;

  89.     /** Time steps from Butcher array (without the first zero). */
  90.     private double[] realC = new double[0];

  92.     /** Internal weights from Butcher array (without the first empty row). Real version, optional. */
  93.     private double[][] realA = new double[0][];

  95.     /** External weights for the high order method from Butcher array. Real version, optional. */
  96.     private double[] realB = new double[0];

  98.     /** Stepsize control exponent. */
  99.     private final double exp;

  101.     /** Safety factor for stepsize control. */
  102.     private T safety;

  104.     /** Minimal reduction factor for stepsize control. */
  105.     private T minReduction;

  107.     /** Maximal growth factor for stepsize control. */
  108.     private T maxGrowth;

  110.     /** Flag setting whether coefficients in Butcher array are interpreted as Field or real numbers. */
  111.     private boolean usingFieldCoefficients;

  113.     /** Build a Runge-Kutta integrator with the given Butcher array.
  114.      * @param field field to which the time and state vector elements belong
  115.      * @param name name of the method
  116.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  117.      * or -1 if method is not <i>fsal</i>
  118.      * @param minStep minimal step (sign is irrelevant, regardless of
  119.      * integration direction, forward or backward), the last step can
  120.      * be smaller than this
  121.      * @param maxStep maximal step (sign is irrelevant, regardless of
  122.      * integration direction, forward or backward), the last step can
  123.      * be smaller than this
  124.      * @param scalAbsoluteTolerance allowed absolute error
  125.      * @param scalRelativeTolerance allowed relative error
  126.      */
  127.     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
  128.                                                 final double minStep, final double maxStep,
  129.                                                 final double scalAbsoluteTolerance,
  130.                                                 final double scalRelativeTolerance) {

  132.         super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

  134.         this.fsal = fsal;
  135.         this.c    = getC();
  136.         this.a    = getA();
  137.         this.b    = getB();

  139.         exp = -1.0 / getOrder();

  141.         // set the default values of the algorithm control parameters
  142.         setSafety(field.getZero().add(0.9));
  143.         setMinReduction(field.getZero().add(0.2));
  144.         setMaxGrowth(field.getZero().add(10.0));

  146.     }

  148.     /** Build a Runge-Kutta integrator with the given Butcher array.
  149.      * @param field field to which the time and state vector elements belong
  150.      * @param name name of the method
  151.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  152.      * or -1 if method is not <i>fsal</i>
  153.      * @param minStep minimal step (must be positive even for backward
  154.      * integration), the last step can be smaller than this
  155.      * @param maxStep maximal step (must be positive even for backward
  156.      * integration)
  157.      * @param vecAbsoluteTolerance allowed absolute error
  158.      * @param vecRelativeTolerance allowed relative error
  159.      */
  160.     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
  161.                                                 final double   minStep, final double maxStep,
  162.                                                 final double[] vecAbsoluteTolerance,
  163.                                                 final double[] vecRelativeTolerance) {

  165.         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  166.         this.usingFieldCoefficients = false;

  168.         this.fsal = fsal;
  169.         this.c    = getC();
  170.         this.a    = getA();
  171.         this.b    = getB();

  173.         exp = -1.0 / getOrder();

  175.         // set the default values of the algorithm control parameters
  176.         setSafety(field.getZero().add(0.9));
  177.         setMinReduction(field.getZero().add(0.2));
  178.         setMaxGrowth(field.getZero().add(10.0));

  180.     }

  182.     /** Create an interpolator.
  183.      * @param forward integration direction indicator
  184.      * @param yDotK slopes at the intermediate points
  185.      * @param globalPreviousState start of the global step
  186.      * @param globalCurrentState end of the global step
  187.      * @param mapper equations mapper for the all equations
  188.      * @return external weights for the high order method from Butcher array
  189.      */
  190.     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
  191.                                                                               FieldODEStateAndDerivative<T> globalPreviousState,
  192.                                                                               FieldODEStateAndDerivative<T> globalCurrentState,
  193.                                                                               FieldEquationsMapper<T> mapper);

  195.     /** Get the order of the method.
  196.      * @return order of the method
  197.      */
  198.     public abstract int getOrder();

  200.     /** Get the safety factor for stepsize control.
  201.      * @return safety factor
  202.      */
  203.     public T getSafety() {
  204.         return safety;
  205.     }

  207.     /** Set the safety factor for stepsize control.
  208.      * @param safety safety factor
  209.      */
  210.     public void setSafety(final T safety) {
  211.         this.safety = safety;
  212.     }

  214.     /**
  215.      * Setter for the flag between real or Field coefficients in the Butcher array.
  216.      *
  217.      * @param usingFieldCoefficients new value for flag
  218.      */
  219.     public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
  220.         this.usingFieldCoefficients = usingFieldCoefficients;
  221.     }

  223.     /** {@inheritDoc} */
  224.     @Override
  225.     public boolean isUsingFieldCoefficients() {
  226.         return usingFieldCoefficients;
  227.     }

  229.     /** {@inheritDoc} */
  230.     @Override
  231.     public int getNumberOfStages() {
  232.         return b.length;
  233.     }

  235.     /** {@inheritDoc} */
  236.     @Override
  237.     protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
  238.         if (!isUsingFieldCoefficients()) {
  239.             realA = getRealA();
  240.             realB = getRealB();
  241.             realC = getRealC();
  242.         }
  243.         return super.initIntegration(eqn, s0, t);
  244.     }

  246.     /** {@inheritDoc} */
  247.     @Override
  248.     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
  249.                                                    final FieldODEState<T> initialState, final T finalTime)
  250.         throws MathIllegalArgumentException, MathIllegalStateException {

  252.         sanityChecks(initialState, finalTime);
  253.         setStepStart(initIntegration(equations, initialState, finalTime));
  254.         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;

  256.         // create some internal working arrays
  257.         final int   stages = getNumberOfStages();
  258.         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
  259.         T[]   yTmp   = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());

  261.         // set up integration control objects
  262.         T  hNew           = getField().getZero();
  263.         boolean firstTime = true;

  265.         // main integration loop
  266.         setIsLastStep(false);
  267.         do {

  269.             // iterate over step size, ensuring local normalized error is smaller than 1
  270.             double error = 10.0;
  271.             while (error >= 1.0) {

  273.                 // first stage
  274.                 final T[] y = getStepStart().getCompleteState();
  275.                 yDotK[0] = getStepStart().getCompleteDerivative();

  277.                 if (firstTime) {
  278.                     final StepsizeHelper helper = getStepSizeHelper();
  279.                     final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
  280.                     for (int i = 0; i < scale.length; ++i) {
  281.                         scale[i] = helper.getTolerance(i, y[i].abs());
  282.                     }
  283.                     hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart()));
  284.                     firstTime = false;
  285.                 }

  287.                 setStepSize(hNew);
  288.                 if (forward) {
  289.                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
  290.                         setStepSize(finalTime.subtract(getStepStart().getTime()));
  291.                     }
  292.                 } else {
  293.                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
  294.                         setStepSize(finalTime.subtract(getStepStart().getTime()));
  295.                     }
  296.                 }

  298.                 // next stages
  299.                 if (isUsingFieldCoefficients()) {
  300.                     FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
  301.                             getStepStart().getTime(), y, getStepSize(), a, c, yDotK);
  302.                     yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
  303.                 } else {
  304.                     FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
  305.                             getStepStart().getTime(), y, getStepSize(), realA, realC, yDotK);
  306.                     yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
  307.                 }

  309.                 incrementEvaluations(stages - 1);

  311.                 // estimate the error at the end of the step
  312.                 error = estimateError(yDotK, y, yTmp, getStepSize());
  313.                 if (error >= 1.0) {
  314.                     // reject the step and attempt to reduce error by stepsize control
  315.                     final T factor = MathUtils.min(maxGrowth,
  316.                                                    MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
  317.                     hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
  318.                 }

  320.             }
  321.             final T   stepEnd = getStepStart().getTime().add(getStepSize());
  322.             final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
  323.             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);

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

  329.             if (!isLastStep()) {

  331.                 // stepsize control for next step
  332.                 final T factor = MathUtils.min(maxGrowth,
  333.                                                MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
  334.                 final T  scaledH    = getStepSize().multiply(factor);
  335.                 final T  nextT      = getStepStart().getTime().add(scaledH);
  336.                 final boolean nextIsLast = forward ?
  337.                                            nextT.subtract(finalTime).getReal() >= 0 :
  338.                                            nextT.subtract(finalTime).getReal() <= 0;
  339.                 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);

  341.                 final T  filteredNextT      = getStepStart().getTime().add(hNew);
  342.                 final boolean filteredNextIsLast = forward ?
  343.                                                    filteredNextT.subtract(finalTime).getReal() >= 0 :
  344.                                                    filteredNextT.subtract(finalTime).getReal() <= 0;
  345.                 if (filteredNextIsLast) {
  346.                     hNew = finalTime.subtract(getStepStart().getTime());
  347.                 }

  349.             }

  351.         } while (!isLastStep());

  353.         final FieldODEStateAndDerivative<T> finalState = getStepStart();
  354.         resetInternalState();
  355.         return finalState;

  357.     }

  359.     /** Get the minimal reduction factor for stepsize control.
  360.      * @return minimal reduction factor
  361.      */
  362.     public T getMinReduction() {
  363.         return minReduction;
  364.     }

  366.     /** Set the minimal reduction factor for stepsize control.
  367.      * @param minReduction minimal reduction factor
  368.      */
  369.     public void setMinReduction(final T minReduction) {
  370.         this.minReduction = minReduction;
  371.     }

  373.     /** Get the maximal growth factor for stepsize control.
  374.      * @return maximal growth factor
  375.      */
  376.     public T getMaxGrowth() {
  377.         return maxGrowth;
  378.     }

  380.     /** Set the maximal growth factor for stepsize control.
  381.      * @param maxGrowth maximal growth factor
  382.      */
  383.     public void setMaxGrowth(final T maxGrowth) {
  384.         this.maxGrowth = maxGrowth;
  385.     }

  387.     /** Compute the error ratio.
  388.      * @param yDotK derivatives computed during the first stages
  389.      * @param y0 estimate of the step at the start of the step
  390.      * @param y1 estimate of the step at the end of the step
  391.      * @param h  current step
  392.      * @return error ratio, greater than 1 if step should be rejected
  393.      */
  394.     protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);

  396. }
Created with JaCoCo 0.8.12.202403310830