MultistepFieldIntegrator.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.Field;
  26. import org.hipparchus.CalculusFieldElement;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.exception.MathIllegalStateException;
  29. import org.hipparchus.linear.Array2DRowFieldMatrix;
  30. import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
  31. import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
  32. import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
  33. import org.hipparchus.ode.sampling.FieldODEStepHandler;
  34. import org.hipparchus.util.FastMath;
  35. import org.hipparchus.util.MathArrays;

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

  72.     /** First scaled derivative (h y'). */
  73.     protected T[] scaled;

  75.     /** Nordsieck matrix of the higher scaled derivatives.
  76.      * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
  77.      */
  78.     protected Array2DRowFieldMatrix<T> nordsieck;

  80.     /** Starter integrator. */
  81.     private FieldODEIntegrator<T> starter;

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

  86.     /** Stepsize control exponent. */
  87.     private double exp;

  89.     /** Safety factor for stepsize control. */
  90.     private double safety;

  92.     /** Minimal reduction factor for stepsize control. */
  93.     private double minReduction;

  95.     /** Maximal growth factor for stepsize control. */
  96.     private double maxGrowth;

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

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

  128.         if (nSteps < 2) {
  129.             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  130.                                                    nSteps, 2, true);
  131.         }

  133.         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
  134.                                                         scalAbsoluteTolerance,
  135.                                                         scalRelativeTolerance);
  136.         this.nSteps = nSteps;

  138.         exp = -1.0 / order;

  140.         // set the default values of the algorithm control parameters
  141.         setSafety(0.9);
  142.         setMinReduction(0.2);
  143.         setMaxGrowth(FastMath.pow(2.0, -exp));

  145.     }

  147.     /**
  148.      * Build a multistep integrator with the given stepsize bounds.
  149.      * <p>The default starter integrator is set to the {@link
  150.      * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
  151.      * some defaults settings.</p>
  152.      * <p>
  153.      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
  154.      * </p>
  155.      * @param field field to which the time and state vector elements belong
  156.      * @param name name of the method
  157.      * @param nSteps number of steps of the multistep method
  158.      * (excluding the one being computed)
  159.      * @param order order of the method
  160.      * @param minStep minimal step (must be positive even for backward
  161.      * integration), the last step can be smaller than this
  162.      * @param maxStep maximal step (must be positive even for backward
  163.      * integration)
  164.      * @param vecAbsoluteTolerance allowed absolute error
  165.      * @param vecRelativeTolerance allowed relative error
  166.      */
  167.     protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
  168.                                        final int order,
  169.                                        final double minStep, final double maxStep,
  170.                                        final double[] vecAbsoluteTolerance,
  171.                                        final double[] vecRelativeTolerance) {
  172.         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);

  174.         if (nSteps < 2) {
  175.             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  176.                                                    nSteps, 2, true);
  177.         }

  179.         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
  180.                                                         vecAbsoluteTolerance,
  181.                                                         vecRelativeTolerance);
  182.         this.nSteps = nSteps;

  184.         exp = -1.0 / order;

  186.         // set the default values of the algorithm control parameters
  187.         setSafety(0.9);
  188.         setMinReduction(0.2);
  189.         setMaxGrowth(FastMath.pow(2.0, -exp));

  191.     }

  193.     /**
  194.      * Get the starter integrator.
  195.      * @return starter integrator
  196.      */
  197.     public FieldODEIntegrator<T> getStarterIntegrator() {
  198.         return starter;
  199.     }

  201.     /**
  202.      * Set the starter integrator.
  203.      * <p>The various step and event handlers for this starter integrator
  204.      * will be managed automatically by the multi-step integrator. Any
  205.      * user configuration for these elements will be cleared before use.</p>
  206.      * @param starterIntegrator starter integrator
  207.      */
  208.     public void setStarterIntegrator(FieldODEIntegrator<T> starterIntegrator) {
  209.         this.starter = starterIntegrator;
  210.     }

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

  233.         // make sure NO user events nor user step handlers are triggered,
  234.         // this is the task of the top level integrator, not the task of the starter integrator
  235.         starter.clearEventDetectors();
  236.         starter.clearStepHandlers();

  238.         // set up one specific step handler to extract initial Nordsieck vector
  239.         starter.addStepHandler(new FieldNordsieckInitializer((nSteps + 3) / 2));

  241.         // start integration, expecting a InitializationCompletedMarkerException
  242.         try {

  244.             starter.integrate(equations, initialState, t);

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

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

  252.             // count the evaluations used by the starter
  253.             getEvaluationsCounter().increment(starter.getEvaluations());

  255.         }

  257.         // remove the specific step handler
  258.         starter.clearStepHandlers();

  260.     }

  262.     /** Initialize the high order scaled derivatives at step start.
  263.      * @param h step size to use for scaling
  264.      * @param t first steps times
  265.      * @param y first steps states
  266.      * @param yDot first steps derivatives
  267.      * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
  268.      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
  269.      */
  270.     protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(T h, T[] t, T[][] y, T[][] yDot);

  272.     /** Get the minimal reduction factor for stepsize control.
  273.      * @return minimal reduction factor
  274.      */
  275.     public double getMinReduction() {
  276.         return minReduction;
  277.     }

  279.     /** Set the minimal reduction factor for stepsize control.
  280.      * @param minReduction minimal reduction factor
  281.      */
  282.     public void setMinReduction(final double minReduction) {
  283.         this.minReduction = minReduction;
  284.     }

  286.     /** Get the maximal growth factor for stepsize control.
  287.      * @return maximal growth factor
  288.      */
  289.     public double getMaxGrowth() {
  290.         return maxGrowth;
  291.     }

  293.     /** Set the maximal growth factor for stepsize control.
  294.      * @param maxGrowth maximal growth factor
  295.      */
  296.     public void setMaxGrowth(final double maxGrowth) {
  297.         this.maxGrowth = maxGrowth;
  298.     }

  300.     /** Get the safety factor for stepsize control.
  301.      * @return safety factor
  302.      */
  303.     public double getSafety() {
  304.       return safety;
  305.     }

  307.     /** Set the safety factor for stepsize control.
  308.      * @param safety safety factor
  309.      */
  310.     public void setSafety(final double safety) {
  311.       this.safety = safety;
  312.     }

  314.     /** Get the number of steps of the multistep method (excluding the one being computed).
  315.      * @return number of steps of the multistep method (excluding the one being computed)
  316.      */
  317.     public int getNSteps() {
  318.       return nSteps;
  319.     }

  321.     /** Rescale the instance.
  322.      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
  323.      * this method has the side effect of rescaling this arrays in the caller too.</p>
  324.      * @param newStepSize new step size to use in the scaled and Nordsieck arrays
  325.      */
  326.     protected void rescale(final T newStepSize) {

  328.         final T ratio = newStepSize.divide(getStepSize());
  329.         for (int i = 0; i < scaled.length; ++i) {
  330.             scaled[i] = scaled[i].multiply(ratio);
  331.         }

  333.         final T[][] nData = nordsieck.getDataRef();
  334.         T power = ratio;
  335.         for (int i = 0; i < nData.length; ++i) {
  336.             power = power.multiply(ratio);
  337.             final T[] nDataI = nData[i];
  338.             for (int j = 0; j < nDataI.length; ++j) {
  339.                 nDataI[j] = nDataI[j].multiply(power);
  340.             }
  341.         }

  343.         setStepSize(newStepSize);

  345.     }

  347.     /** Compute step grow/shrink factor according to normalized error.
  348.      * @param error normalized error of the current step
  349.      * @return grow/shrink factor for next step
  350.      */
  351.     protected double computeStepGrowShrinkFactor(final double error) {
  352.         return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
  353.     }

  355.     /** Specialized step handler storing the first step.
  356.      */
  357.     private class FieldNordsieckInitializer implements FieldODEStepHandler<T> {

  359.         /** Steps counter. */
  360.         private int count;

  362.         /** Saved start. */
  363.         private FieldODEStateAndDerivative<T> savedStart;

  365.         /** First steps times. */
  366.         private final T[] t;

  368.         /** First steps states. */
  369.         private final T[][] y;

  371.         /** First steps derivatives. */
  372.         private final T[][] yDot;

  374.         /** Simple constructor.
  375.          * @param nbStartPoints number of start points (including the initial point)
  376.          */
  377.         FieldNordsieckInitializer(final int nbStartPoints) {
  378.             this.count  = 0;
  379.             this.t      = MathArrays.buildArray(getField(), nbStartPoints);
  380.             this.y      = MathArrays.buildArray(getField(), nbStartPoints, -1);
  381.             this.yDot   = MathArrays.buildArray(getField(), nbStartPoints, -1);
  382.         }

  384.         /** {@inheritDoc} */
  385.         @Override
  386.         public void handleStep(FieldODEStateInterpolator<T> interpolator) {


  389.             if (count == 0) {
  390.                 // first step, we need to store also the point at the beginning of the step
  391.                 final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
  392.                 savedStart  = prev;
  393.                 t[count]    = prev.getTime();
  394.                 y[count]    = prev.getCompleteState();
  395.                 yDot[count] = prev.getCompleteDerivative();
  396.             }

  398.             // store the point at the end of the step
  399.             ++count;
  400.             final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
  401.             t[count]    = curr.getTime();
  402.             y[count]    = curr.getCompleteState();
  403.             yDot[count] = curr.getCompleteDerivative();

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

  407.                 // this was the last point we needed, we can compute the derivatives
  408.                 setStepStart(savedStart);
  409.                 final T rawStep = t[t.length - 1].subtract(t[0]).divide(t.length - 1);
  410.                 setStepSize(getStepSizeHelper().filterStep(rawStep, rawStep.getReal() >= 0, true));

  412.                 // first scaled derivative
  413.                 scaled = MathArrays.buildArray(getField(), yDot[0].length);
  414.                 for (int j = 0; j < scaled.length; ++j) {
  415.                     scaled[j] = yDot[0][j].multiply(getStepSize());
  416.                 }

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

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

  424.             }

  426.         }

  428.         /** {@inheritDoc} */
  429.         @Override
  430.         public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
  431.             // nothing to do
  432.         }

  434.     }

  436.     /** Marker exception used ONLY to stop the starter integrator after first step. */
  437.     private static class InitializationCompletedMarkerException
  438.         extends RuntimeException {

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

  443.         /** Simple constructor. */
  444.         InitializationCompletedMarkerException() {
  445.             super((Throwable) null);
  446.         }

  448.     }

  450. }
Created with JaCoCo 0.8.11.202310140853