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

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

  20. import org.apache.commons.math4.legacy.core.Field;
  21. import org.apache.commons.math4.legacy.core.RealFieldElement;
  22. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  23. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  24. import org.apache.commons.math4.legacy.exception.NoBracketingException;
  25. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  26. import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
  27. import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
  28. import org.apache.commons.math4.legacy.ode.FieldODEState;
  29. import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
  30. import org.apache.commons.math4.legacy.core.MathArrays;

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

  68. public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
  69.     extends AdaptiveStepsizeFieldIntegrator<T>
  70.     implements FieldButcherArrayProvider<T> {

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

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

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

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

  84.     /** Stepsize control exponent. */
  85.     private final T exp;

  87.     /** Safety factor for stepsize control. */
  88.     private T safety;

  90.     /** Minimal reduction factor for stepsize control. */
  91.     private T minReduction;

  93.     /** Maximal growth factor for stepsize control. */
  94.     private T maxGrowth;

  96.     /** Build a Runge-Kutta integrator with the given Butcher array.
  97.      * @param field field to which the time and state vector elements belong
  98.      * @param name name of the method
  99.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  100.      * or -1 if method is not <i>fsal</i>
  101.      * @param minStep minimal step (sign is irrelevant, regardless of
  102.      * integration direction, forward or backward), the last step can
  103.      * be smaller than this
  104.      * @param maxStep maximal step (sign is irrelevant, regardless of
  105.      * integration direction, forward or backward), the last step can
  106.      * be smaller than this
  107.      * @param scalAbsoluteTolerance allowed absolute error
  108.      * @param scalRelativeTolerance allowed relative error
  109.      */
  110.     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
  111.                                                 final double minStep, final double maxStep,
  112.                                                 final double scalAbsoluteTolerance,
  113.                                                 final double scalRelativeTolerance) {

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

  117.         this.fsal = fsal;
  118.         this.c    = getC();
  119.         this.a    = getA();
  120.         this.b    = getB();

  122.         exp = field.getOne().divide(-getOrder());

  124.         // set the default values of the algorithm control parameters
  125.         setSafety(field.getZero().add(0.9));
  126.         setMinReduction(field.getZero().add(0.2));
  127.         setMaxGrowth(field.getZero().add(10.0));
  128.     }

  130.     /** Build a Runge-Kutta integrator with the given Butcher array.
  131.      * @param field field to which the time and state vector elements belong
  132.      * @param name name of the method
  133.      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
  134.      * or -1 if method is not <i>fsal</i>
  135.      * @param minStep minimal step (must be positive even for backward
  136.      * integration), the last step can be smaller than this
  137.      * @param maxStep maximal step (must be positive even for backward
  138.      * integration)
  139.      * @param vecAbsoluteTolerance allowed absolute error
  140.      * @param vecRelativeTolerance allowed relative error
  141.      */
  142.     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
  143.                                                 final double   minStep, final double maxStep,
  144.                                                 final double[] vecAbsoluteTolerance,
  145.                                                 final double[] vecRelativeTolerance) {

  147.         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);

  149.         this.fsal = fsal;
  150.         this.c    = getC();
  151.         this.a    = getA();
  152.         this.b    = getB();

  154.         exp = field.getOne().divide(-getOrder());

  156.         // set the default values of the algorithm control parameters
  157.         setSafety(field.getZero().add(0.9));
  158.         setMinReduction(field.getZero().add(0.2));
  159.         setMaxGrowth(field.getZero().add(10.0));
  160.     }

  162.     /** Create a fraction.
  163.      * @param p numerator
  164.      * @param q denominator
  165.      * @return p/q computed in the instance field
  166.      */
  167.     protected T fraction(final int p, final int q) {
  168.         return getField().getOne().multiply(p).divide(q);
  169.     }

  171.     /** Create a fraction.
  172.      * @param p numerator
  173.      * @param q denominator
  174.      * @return p/q computed in the instance field
  175.      */
  176.     protected T fraction(final double p, final double q) {
  177.         return getField().getOne().multiply(p).divide(q);
  178.     }

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

  197.     /** Get the safety factor for stepsize control.
  198.      * @return safety factor
  199.      */
  200.     public T getSafety() {
  201.         return safety;
  202.     }

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

  211.     /** {@inheritDoc} */
  212.     @Override
  213.     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
  214.                                                    final FieldODEState<T> initialState, final T finalTime)
  215.         throws NumberIsTooSmallException, DimensionMismatchException,
  216.         MaxCountExceededException, NoBracketingException {

  218.         sanityChecks(initialState, finalTime);
  219.         final T   t0 = initialState.getTime();
  220.         final T[] y0 = equations.getMapper().mapState(initialState);
  221.         setStepStart(initIntegration(equations, t0, y0, finalTime));
  222.         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;

  224.         // create some internal working arrays
  225.         final int   stages = c.length + 1;
  226.         T[]         y      = y0;
  227.         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
  228.         final T[]   yTmp   = MathArrays.buildArray(getField(), y0.length);

  230.         // set up integration control objects
  231.         T  hNew           = getField().getZero();
  232.         boolean firstTime = true;

  234.         // main integration loop
  235.         setIsLastStep(false);
  236.         do {

  238.             // iterate over step size, ensuring local normalized error is smaller than 1
  239.             T error = getField().getZero().add(10);
  240.             while (error.subtract(1.0).getReal() >= 0) {

  242.                 // first stage
  243.                 y        = equations.getMapper().mapState(getStepStart());
  244.                 yDotK[0] = equations.getMapper().mapDerivative(getStepStart());

  246.                 if (firstTime) {
  247.                     final T[] scale = MathArrays.buildArray(getField(), mainSetDimension);
  248.                     if (vecAbsoluteTolerance == null) {
  249.                         for (int i = 0; i < scale.length; ++i) {
  250.                             scale[i] = y[i].abs().multiply(scalRelativeTolerance).add(scalAbsoluteTolerance);
  251.                         }
  252.                     } else {
  253.                         for (int i = 0; i < scale.length; ++i) {
  254.                             scale[i] = y[i].abs().multiply(vecRelativeTolerance[i]).add(vecAbsoluteTolerance[i]);
  255.                         }
  256.                     }
  257.                     hNew = initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper());
  258.                     firstTime = false;
  259.                 }

  261.                 setStepSize(hNew);
  262.                 if (forward) {
  263.                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
  264.                         setStepSize(finalTime.subtract(getStepStart().getTime()));
  265.                     }
  266.                 } else {
  267.                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
  268.                         setStepSize(finalTime.subtract(getStepStart().getTime()));
  269.                     }
  270.                 }

  272.                 // next stages
  273.                 for (int k = 1; k < stages; ++k) {

  275.                     for (int j = 0; j < y0.length; ++j) {
  276.                         T sum = yDotK[0][j].multiply(a[k-1][0]);
  277.                         for (int l = 1; l < k; ++l) {
  278.                             sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
  279.                         }
  280.                         yTmp[j] = y[j].add(getStepSize().multiply(sum));
  281.                     }

  283.                     yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
  284.                 }

  286.                 // estimate the state at the end of the step
  287.                 for (int j = 0; j < y0.length; ++j) {
  288.                     T sum    = yDotK[0][j].multiply(b[0]);
  289.                     for (int l = 1; l < stages; ++l) {
  290.                         sum = sum.add(yDotK[l][j].multiply(b[l]));
  291.                     }
  292.                     yTmp[j] = y[j].add(getStepSize().multiply(sum));
  293.                 }

  295.                 // estimate the error at the end of the step
  296.                 error = estimateError(yDotK, y, yTmp, getStepSize());
  297.                 if (error.subtract(1.0).getReal() >= 0) {
  298.                     // reject the step and attempt to reduce error by stepsize control
  299.                     final T factor = RealFieldElement.min(maxGrowth,
  300.                                                    RealFieldElement.max(minReduction, safety.multiply(error.pow(exp))));
  301.                     hNew = filterStep(getStepSize().multiply(factor), forward, false);
  302.                 }
  303.             }
  304.             final T   stepEnd = getStepStart().getTime().add(getStepSize());
  305.             final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
  306.             final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<>(stepEnd, yTmp, yDotTmp);

  308.             // local error is small enough: accept the step, trigger events and step handlers
  309.             System.arraycopy(yTmp, 0, y, 0, y0.length);
  310.             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
  311.                                     finalTime));

  313.             if (!isLastStep()) {

  315.                 // stepsize control for next step
  316.                 final T factor = RealFieldElement.min(maxGrowth,
  317.                                                RealFieldElement.max(minReduction, safety.multiply(error.pow(exp))));
  318.                 final T  scaledH    = getStepSize().multiply(factor);
  319.                 final T  nextT      = getStepStart().getTime().add(scaledH);
  320.                 final boolean nextIsLast = forward ?
  321.                                            nextT.subtract(finalTime).getReal() >= 0 :
  322.                                            nextT.subtract(finalTime).getReal() <= 0;
  323.                 hNew = filterStep(scaledH, forward, nextIsLast);

  325.                 final T  filteredNextT      = getStepStart().getTime().add(hNew);
  326.                 final boolean filteredNextIsLast = forward ?
  327.                                                    filteredNextT.subtract(finalTime).getReal() >= 0 :
  328.                                                    filteredNextT.subtract(finalTime).getReal() <= 0;
  329.                 if (filteredNextIsLast) {
  330.                     hNew = finalTime.subtract(getStepStart().getTime());
  331.                 }
  332.             }
  333.         } while (!isLastStep());

  335.         final FieldODEStateAndDerivative<T> finalState = getStepStart();
  336.         resetInternalState();
  337.         return finalState;
  338.     }

  340.     /** Get the minimal reduction factor for stepsize control.
  341.      * @return minimal reduction factor
  342.      */
  343.     public T getMinReduction() {
  344.         return minReduction;
  345.     }

  347.     /** Set the minimal reduction factor for stepsize control.
  348.      * @param minReduction minimal reduction factor
  349.      */
  350.     public void setMinReduction(final T minReduction) {
  351.         this.minReduction = minReduction;
  352.     }

  354.     /** Get the maximal growth factor for stepsize control.
  355.      * @return maximal growth factor
  356.      */
  357.     public T getMaxGrowth() {
  358.         return maxGrowth;
  359.     }

  361.     /** Set the maximal growth factor for stepsize control.
  362.      * @param maxGrowth maximal growth factor
  363.      */
  364.     public void setMaxGrowth(final T maxGrowth) {
  365.         this.maxGrowth = maxGrowth;
  366.     }

  368.     /** Compute the error ratio.
  369.      * @param yDotK derivatives computed during the first stages
  370.      * @param y0 estimate of the step at the start of the step
  371.      * @param y1 estimate of the step at the end of the step
  372.      * @param h  current step
  373.      * @return error ratio, greater than 1 if step should be rejected
  374.      */
  375.     protected abstract T estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
  376. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT