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

  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathIllegalStateException;
  27. import org.hipparchus.linear.Array2DRowRealMatrix;
  28. import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
  29. import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
  30. import org.hipparchus.ode.sampling.ODEStateInterpolator;
  31. import org.hipparchus.ode.sampling.ODEStepHandler;
  32. import org.hipparchus.util.FastMath;

  34. /**
  35.  * This class is the base class for multistep integrators for Ordinary
  36.  * Differential Equations.
  37.  * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
  38.  * \[
  39.  *   \left\{\begin{align}
  40.  *   s_1(n) &amp;= h y'_n \text{ for first derivative}\\
  41.  *   s_2(n) &amp;= \frac{h^2}{2} y_n'' \text{ for second derivative}\\
  42.  *   s_3(n) &amp;= \frac{h^3}{6} y_n''' \text{ for third derivative}\\
  43.  *   &amp;\cdots\\
  44.  *   s_k(n) &amp;= \frac{h^k}{k!} y_n^{(k)} \text{ for } k^\mathrm{th} \text{ derivative}
  45.  *   \end{align}\right.
  46.  * \]</p>
  47.  * <p>Rather than storing several previous steps separately, this implementation uses
  48.  * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
  49.  * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
  50.  * \[
  51.  * r_n = [ s_2(n), s_3(n) \ldots s_k(n) ]^T
  52.  * \]
  53.  * (we omit the k index in the notation for clarity)</p>
  54.  * <p>
  55.  * Multistep integrators with Nordsieck representation are highly sensitive to
  56.  * large step changes because when the step is multiplied by factor a, the
  57.  * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
  58.  * and the last components are the least accurate ones. The default max growth
  59.  * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
  60.  * </p>
  61.  *
  62.  * @see org.hipparchus.ode.nonstiff.AdamsBashforthIntegrator
  63.  * @see org.hipparchus.ode.nonstiff.AdamsMoultonIntegrator
  64.  */
  65. public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {

  67.     /** First scaled derivative (h y'). */
  68.     protected double[] scaled;

  70.     /** Nordsieck matrix of the higher scaled derivatives.
  71.      * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
  72.      */
  73.     protected Array2DRowRealMatrix nordsieck;

  75.     /** Starter integrator. */
  76.     private ODEIntegrator starter;

  78.     /** Number of steps of the multistep method (excluding the one being computed). */
  79.     private final int nSteps;

  81.     /** Stepsize control exponent. */
  82.     private double exp;

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

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

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

  93.     /**
  94.      * Build a multistep integrator with the given stepsize bounds.
  95.      * <p>The default starter integrator is set to the {@link
  96.      * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
  97.      * some defaults settings.</p>
  98.      * <p>
  99.      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
  100.      * </p>
  101.      * @param name name of the method
  102.      * @param nSteps number of steps of the multistep method
  103.      * (excluding the one being computed)
  104.      * @param order order of the method
  105.      * @param minStep minimal step (must be positive even for backward
  106.      * integration), the last step can be smaller than this
  107.      * @param maxStep maximal step (must be positive even for backward
  108.      * integration)
  109.      * @param scalAbsoluteTolerance allowed absolute error
  110.      * @param scalRelativeTolerance allowed relative error
  111.      * @exception MathIllegalArgumentException if number of steps is smaller than 2
  112.      */
  113.     protected MultistepIntegrator(final String name, final int nSteps,
  114.                                   final int order,
  115.                                   final double minStep, final double maxStep,
  116.                                   final double scalAbsoluteTolerance,
  117.                                   final double scalRelativeTolerance)
  118.         throws MathIllegalArgumentException {

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

  122.         if (nSteps < 2) {
  123.             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  124.                                                    nSteps, 2, true);
  125.         }

  127.         starter = new DormandPrince853Integrator(minStep, maxStep,
  128.                                                  scalAbsoluteTolerance,
  129.                                                  scalRelativeTolerance);
  130.         this.nSteps = nSteps;

  132.         exp = -1.0 / order;

  134.         // set the default values of the algorithm control parameters
  135.         setSafety(0.9);
  136.         setMinReduction(0.2);
  137.         setMaxGrowth(FastMath.pow(2.0, -exp));

  139.     }

  141.     /**
  142.      * Build a multistep integrator with the given stepsize bounds.
  143.      * <p>The default starter integrator is set to the {@link
  144.      * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
  145.      * some defaults settings.</p>
  146.      * <p>
  147.      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
  148.      * </p>
  149.      * @param name name of the method
  150.      * @param nSteps number of steps of the multistep method
  151.      * (excluding the one being computed)
  152.      * @param order order of the method
  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 MultistepIntegrator(final String name, final int nSteps,
  161.                                   final int order,
  162.                                   final double minStep, final double maxStep,
  163.                                   final double[] vecAbsoluteTolerance,
  164.                                   final double[] vecRelativeTolerance) {
  165.         super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);

  167.         if (nSteps < 2) {
  168.             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  169.                                                    nSteps, 2, true);
  170.         }

  172.         starter = new DormandPrince853Integrator(minStep, maxStep,
  173.                                                  vecAbsoluteTolerance,
  174.                                                  vecRelativeTolerance);
  175.         this.nSteps = nSteps;

  177.         exp = -1.0 / order;

  179.         // set the default values of the algorithm control parameters
  180.         setSafety(0.9);
  181.         setMinReduction(0.2);
  182.         setMaxGrowth(FastMath.pow(2.0, -exp));

  184.     }

  186.     /**
  187.      * Get the starter integrator.
  188.      * @return starter integrator
  189.      */
  190.     public ODEIntegrator getStarterIntegrator() {
  191.         return starter;
  192.     }

  194.     /**
  195.      * Set the starter integrator.
  196.      * <p>The various step and event handlers for this starter integrator
  197.      * will be managed automatically by the multi-step integrator. Any
  198.      * user configuration for these elements will be cleared before use.</p>
  199.      * @param starterIntegrator starter integrator
  200.      */
  201.     public void setStarterIntegrator(ODEIntegrator starterIntegrator) {
  202.         this.starter = starterIntegrator;
  203.     }

  205.     /** Start the integration.
  206.      * <p>This method computes one step using the underlying starter integrator,
  207.      * and initializes the Nordsieck vector at step start. The starter integrator
  208.      * purpose is only to establish initial conditions, it does not really change
  209.      * time by itself. The top level multistep integrator remains in charge of
  210.      * handling time propagation and events handling as it will starts its own
  211.      * computation right from the beginning. In a sense, the starter integrator
  212.      * can be seen as a dummy one and so it will never trigger any user event nor
  213.      * call any user step handler.</p>
  214.      * @param equations complete set of differential equations to integrate
  215.      * @param initialState initial state (time, primary and secondary state vectors)
  216.      * @param finalTime target time for the integration
  217.      * (can be set to a value smaller than {@code initialState.getTime()} for backward integration)
  218.      * @exception MathIllegalArgumentException if arrays dimension do not match equations settings
  219.      * @exception MathIllegalArgumentException if integration step is too small
  220.      * @exception MathIllegalStateException if the number of functions evaluations is exceeded
  221.      * @exception MathIllegalArgumentException if the location of an event cannot be bracketed
  222.      */
  223.     protected void start(final ExpandableODE equations, final ODEState initialState, final double finalTime)
  224.         throws MathIllegalArgumentException, MathIllegalStateException {

  226.         // make sure NO user events nor user step handlers are triggered,
  227.         // this is the task of the top level integrator, not the task of the starter integrator
  228.         starter.clearEventDetectors();
  229.         starter.clearStepHandlers();

  231.         // set up one specific step handler to extract initial Nordsieck vector
  232.         starter.addStepHandler(new NordsieckInitializer((nSteps + 3) / 2));

  234.         // start integration, expecting a InitializationCompletedMarkerException
  235.         try {

  237.             starter.integrate(getEquations(), initialState, finalTime);

  239.             // we should not reach this step
  240.             throw new MathIllegalStateException(LocalizedODEFormats.MULTISTEP_STARTER_STOPPED_EARLY);

  242.         } catch (InitializationCompletedMarkerException icme) { // NOPMD
  243.             // this is the expected nominal interruption of the start integrator

  245.             // count the evaluations used by the starter
  246.             getEvaluationsCounter().increment(starter.getEvaluations());

  248.         }

  250.         // remove the specific step handler
  251.         starter.clearStepHandlers();

  253.     }

  255.     /** Initialize the high order scaled derivatives at step start.
  256.      * @param h step size to use for scaling
  257.      * @param t first steps times
  258.      * @param y first steps states
  259.      * @param yDot first steps derivatives
  260.      * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
  261.      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
  262.      */
  263.     protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(double h, double[] t, double[][] y, double[][] yDot);

  265.     /** Get the minimal reduction factor for stepsize control.
  266.      * @return minimal reduction factor
  267.      */
  268.     public double getMinReduction() {
  269.         return minReduction;
  270.     }

  272.     /** Set the minimal reduction factor for stepsize control.
  273.      * @param minReduction minimal reduction factor
  274.      */
  275.     public void setMinReduction(final double minReduction) {
  276.         this.minReduction = minReduction;
  277.     }

  279.     /** Get the maximal growth factor for stepsize control.
  280.      * @return maximal growth factor
  281.      */
  282.     public double getMaxGrowth() {
  283.         return maxGrowth;
  284.     }

  286.     /** Set the maximal growth factor for stepsize control.
  287.      * @param maxGrowth maximal growth factor
  288.      */
  289.     public void setMaxGrowth(final double maxGrowth) {
  290.         this.maxGrowth = maxGrowth;
  291.     }

  293.     /** Get the safety factor for stepsize control.
  294.      * @return safety factor
  295.      */
  296.     public double getSafety() {
  297.       return safety;
  298.     }

  300.     /** Set the safety factor for stepsize control.
  301.      * @param safety safety factor
  302.      */
  303.     public void setSafety(final double safety) {
  304.       this.safety = safety;
  305.     }

  307.     /** Get the number of steps of the multistep method (excluding the one being computed).
  308.      * @return number of steps of the multistep method (excluding the one being computed)
  309.      */
  310.     public int getNSteps() {
  311.       return nSteps;
  312.     }

  314.     /** Rescale the instance.
  315.      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
  316.      * this method has the side effect of rescaling this arrays in the caller too.</p>
  317.      * @param newStepSize new step size to use in the scaled and Nordsieck arrays
  318.      */
  319.     protected void rescale(final double newStepSize) {

  321.         final double ratio = newStepSize / getStepSize();
  322.         for (int i = 0; i < scaled.length; ++i) {
  323.             scaled[i] = scaled[i] * ratio;
  324.         }

  326.         final double[][] nData = nordsieck.getDataRef();
  327.         double power = ratio;
  328.         for (int i = 0; i < nData.length; ++i) {
  329.             power = power * ratio;
  330.             final double[] nDataI = nData[i];
  331.             for (int j = 0; j < nDataI.length; ++j) {
  332.                 nDataI[j] = nDataI[j] * power;
  333.             }
  334.         }

  336.         setStepSize(newStepSize);

  338.     }

  340.     /** Compute step grow/shrink factor according to normalized error.
  341.      * @param error normalized error of the current step
  342.      * @return grow/shrink factor for next step
  343.      */
  344.     protected double computeStepGrowShrinkFactor(final double error) {
  345.         return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
  346.     }

  348.     /** Specialized step handler storing the first step. */
  349.     private class NordsieckInitializer implements ODEStepHandler {

  351.         /** Steps counter. */
  352.         private int count;

  354.         /** Start of the integration. */
  355.         private ODEStateAndDerivative savedStart;

  357.         /** First steps times. */
  358.         private final double[] t;

  360.         /** First steps states. */
  361.         private final double[][] y;

  363.         /** First steps derivatives. */
  364.         private final double[][] yDot;

  366.         /** Simple constructor.
  367.          * @param nbStartPoints number of start points (including the initial point)
  368.          */
  369.         NordsieckInitializer(final int nbStartPoints) {
  370.             this.count  = 0;
  371.             this.t      = new double[nbStartPoints];
  372.             this.y      = new double[nbStartPoints][];
  373.             this.yDot   = new double[nbStartPoints][];
  374.         }

  376.         /** {@inheritDoc} */
  377.         @Override
  378.         public void handleStep(ODEStateInterpolator interpolator) {

  380.             if (count == 0) {
  381.                 // first step, we need to store also the point at the beginning of the step
  382.                 savedStart   = interpolator.getPreviousState();
  383.                 t[0]    = savedStart.getTime();
  384.                 y[0]    = savedStart.getCompleteState();
  385.                 yDot[0] = savedStart.getCompleteDerivative();
  386.             }

  388.             // store the point at the end of the step
  389.             ++count;
  390.             final ODEStateAndDerivative curr = interpolator.getCurrentState();
  391.             t[count]    = curr.getTime();
  392.             y[count]    = curr.getCompleteState();
  393.             yDot[count] = curr.getCompleteDerivative();

  395.             if (count == t.length - 1) {

  397.                 // this was the last point we needed, we can compute the derivatives
  398.                 setStepStart(savedStart);
  399.                 final double rawStep = (t[t.length - 1] - t[0]) / (t.length - 1);
  400.                 setStepSize(getStepSizeHelper().filterStep(rawStep, rawStep >= 0, true));

  402.                 // first scaled derivative
  403.                 scaled = yDot[0].clone();
  404.                 for (int j = 0; j < scaled.length; ++j) {
  405.                     scaled[j] *= getStepSize();
  406.                 }

  408.                 // higher order derivatives
  409.                 nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);

  411.                 // stop the integrator now that all needed steps have been handled
  412.                 throw new InitializationCompletedMarkerException();

  414.             }

  416.         }

  418.     }

  420.     /** Marker exception used ONLY to stop the starter integrator after first step. */
  421.     private static class InitializationCompletedMarkerException
  422.         extends RuntimeException {

  424.         /** Serializable version identifier. */
  425.         private static final long serialVersionUID = -1914085471038046418L;

  427.         /** Simple constructor. */
  428.         InitializationCompletedMarkerException() {
  429.             super((Throwable) null);
  430.         }

  432.     }

  434. }
Created with JaCoCo 0.8.11.202310140853