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

  20. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  21. import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
  22. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  23. import org.apache.commons.math4.legacy.exception.NoBracketingException;
  24. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  25. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  26. import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  27. import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeIntegrator;
  28. import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853Integrator;
  29. import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
  30. import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
  31. import org.apache.commons.math4.core.jdkmath.JdkMath;

  33. /**
  34.  * This class is the base class for multistep integrators for Ordinary
  35.  * Differential Equations.
  36.  * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
  37.  * <div style="white-space: pre"><code>
  38.  * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
  39.  * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
  40.  * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
  41.  * ...
  42.  * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
  43.  * </code></div>
  44.  * <p>Rather than storing several previous steps separately, this implementation uses
  45.  * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
  46.  * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
  47.  * <div style="white-space: pre"><code>
  48.  * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
  49.  * </code></div>
  50.  * (we omit the k index in the notation for clarity)
  51.  * <p>
  52.  * Multistep integrators with Nordsieck representation are highly sensitive to
  53.  * large step changes because when the step is multiplied by factor a, the
  54.  * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
  55.  * and the last components are the least accurate ones. The default max growth
  56.  * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
  57.  * </p>
  58.  *
  59.  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator
  60.  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator
  61.  * @since 2.0
  62.  */
  63. public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {

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

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

  73.     /** Starter integrator. */
  74.     private FirstOrderIntegrator starter;

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

  79.     /** Stepsize control exponent. */
  80.     private double exp;

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

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

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

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

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

  120.         if (nSteps < 2) {
  121.             throw new NumberIsTooSmallException(
  122.                   LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  123.                   nSteps, 2, true);
  124.         }

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

  131.         exp = -1.0 / order;

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

  139.     /**
  140.      * Build a multistep integrator with the given stepsize bounds.
  141.      * <p>The default starter integrator is set to the {@link
  142.      * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
  143.      * some defaults settings.</p>
  144.      * <p>
  145.      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
  146.      * </p>
  147.      * @param name name of the method
  148.      * @param nSteps number of steps of the multistep method
  149.      * (excluding the one being computed)
  150.      * @param order order of the method
  151.      * @param minStep minimal step (must be positive even for backward
  152.      * integration), the last step can be smaller than this
  153.      * @param maxStep maximal step (must be positive even for backward
  154.      * integration)
  155.      * @param vecAbsoluteTolerance allowed absolute error
  156.      * @param vecRelativeTolerance allowed relative error
  157.      */
  158.     protected MultistepIntegrator(final String name, final int nSteps,
  159.                                   final int order,
  160.                                   final double minStep, final double maxStep,
  161.                                   final double[] vecAbsoluteTolerance,
  162.                                   final double[] vecRelativeTolerance) {
  163.         super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
  164.         starter = new DormandPrince853Integrator(minStep, maxStep,
  165.                                                  vecAbsoluteTolerance,
  166.                                                  vecRelativeTolerance);
  167.         this.nSteps = nSteps;

  169.         exp = -1.0 / order;

  171.         // set the default values of the algorithm control parameters
  172.         setSafety(0.9);
  173.         setMinReduction(0.2);
  174.         setMaxGrowth(JdkMath.pow(2.0, -exp));
  175.     }

  177.     /**
  178.      * Get the starter integrator.
  179.      * @return starter integrator
  180.      */
  181.     public ODEIntegrator getStarterIntegrator() {
  182.         return starter;
  183.     }

  185.     /**
  186.      * Set the starter integrator.
  187.      * <p>The various step and event handlers for this starter integrator
  188.      * will be managed automatically by the multi-step integrator. Any
  189.      * user configuration for these elements will be cleared before use.</p>
  190.      * @param starterIntegrator starter integrator
  191.      */
  192.     public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
  193.         this.starter = starterIntegrator;
  194.     }

  196.     /** Start the integration.
  197.      * <p>This method computes one step using the underlying starter integrator,
  198.      * and initializes the Nordsieck vector at step start. The starter integrator
  199.      * purpose is only to establish initial conditions, it does not really change
  200.      * time by itself. The top level multistep integrator remains in charge of
  201.      * handling time propagation and events handling as it will starts its own
  202.      * computation right from the beginning. In a sense, the starter integrator
  203.      * can be seen as a dummy one and so it will never trigger any user event nor
  204.      * call any user step handler.</p>
  205.      * @param t0 initial time
  206.      * @param y0 initial value of the state vector at t0
  207.      * @param t target time for the integration
  208.      * (can be set to a value smaller than <code>t0</code> for backward integration)
  209.      * @exception DimensionMismatchException if arrays dimension do not match equations settings
  210.      * @exception NumberIsTooSmallException if integration step is too small
  211.      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  212.      * @exception NoBracketingException if the location of an event cannot be bracketed
  213.      */
  214.     protected void start(final double t0, final double[] y0, final double t)
  215.         throws DimensionMismatchException, NumberIsTooSmallException,
  216.                MaxCountExceededException, NoBracketingException {

  218.         // make sure NO user event nor user step handler is triggered,
  219.         // this is the task of the top level integrator, not the task
  220.         // of the starter integrator
  221.         starter.clearEventHandlers();
  222.         starter.clearStepHandlers();

  224.         // set up one specific step handler to extract initial Nordsieck vector
  225.         starter.addStepHandler(new NordsieckInitializer((nSteps + 3) / 2, y0.length));

  227.         // start integration, expecting a InitializationCompletedMarkerException
  228.         try {

  230.             if (starter instanceof AbstractIntegrator) {
  231.                 ((AbstractIntegrator) starter).integrate(getExpandable(), t);
  232.             } else {
  233.                 starter.integrate(new FirstOrderDifferentialEquations() {

  235.                     /** {@inheritDoc} */
  236.                     @Override
  237.                     public int getDimension() {
  238.                         return getExpandable().getTotalDimension();
  239.                     }

  241.                     /** {@inheritDoc} */
  242.                     @Override
  243.                     public void computeDerivatives(double t, double[] y, double[] yDot) {
  244.                         getExpandable().computeDerivatives(t, y, yDot);
  245.                     }
  246.                 }, t0, y0, t, new double[y0.length]);
  247.             }

  249.             // we should not reach this step
  250.             throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
  251.         } catch (InitializationCompletedMarkerException icme) { // NOPMD
  252.             // this is the expected nominal interruption of the start integrator

  254.             // count the evaluations used by the starter
  255.             getCounter().increment(starter.getEvaluations());
  256.         }

  258.         // remove the specific step handler
  259.         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 Array2DRowRealMatrix initializeHighOrderDerivatives(double h, double[] t,
  271.                                                                            double[][] y,
  272.                                                                            double[][] yDot);

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

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

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

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

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

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

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

  323.     /** Compute step grow/shrink factor according to normalized error.
  324.      * @param error normalized error of the current step
  325.      * @return grow/shrink factor for next step
  326.      */
  327.     protected double computeStepGrowShrinkFactor(final double error) {
  328.         return JdkMath.min(maxGrowth, JdkMath.max(minReduction, safety * JdkMath.pow(error, exp)));
  329.     }

  331.     /** Transformer used to convert the first step to Nordsieck representation.
  332.      * @deprecated as of 3.6 this unused interface is deprecated
  333.      */
  334.     @Deprecated
  335.     public interface NordsieckTransformer {
  336.         /** Initialize the high order scaled derivatives at step start.
  337.          * @param h step size to use for scaling
  338.          * @param t first steps times
  339.          * @param y first steps states
  340.          * @param yDot first steps derivatives
  341.          * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
  342.          * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
  343.          */
  344.         Array2DRowRealMatrix initializeHighOrderDerivatives(double h, double[] t,
  345.                                                             double[][] y,
  346.                                                             double[][] yDot);
  347.     }

  349.     /** Specialized step handler storing the first step. */
  350.     private final class NordsieckInitializer implements StepHandler {

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

  355.         /** First steps times. */
  356.         private final double[] t;

  358.         /** First steps states. */
  359.         private final double[][] y;

  361.         /** First steps derivatives. */
  362.         private final double[][] yDot;

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

  375.         /** {@inheritDoc} */
  376.         @Override
  377.         public void handleStep(StepInterpolator interpolator, boolean isLast)
  378.             throws MaxCountExceededException {

  380.             final double prev = interpolator.getPreviousTime();
  381.             final double curr = interpolator.getCurrentTime();

  383.             if (count == 0) {
  384.                 // first step, we need to store also the point at the beginning of the step
  385.                 interpolator.setInterpolatedTime(prev);
  386.                 t[0] = prev;
  387.                 final ExpandableStatefulODE expandable = getExpandable();
  388.                 final EquationsMapper primary = expandable.getPrimaryMapper();
  389.                 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
  390.                 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
  391.                 int index = 0;
  392.                 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
  393.                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
  394.                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
  395.                     ++index;
  396.                 }
  397.             }

  399.             // store the point at the end of the step
  400.             ++count;
  401.             interpolator.setInterpolatedTime(curr);
  402.             t[count] = curr;

  404.             final ExpandableStatefulODE expandable = getExpandable();
  405.             final EquationsMapper primary = expandable.getPrimaryMapper();
  406.             primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
  407.             primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
  408.             int index = 0;
  409.             for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
  410.                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
  411.                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
  412.                 ++index;
  413.             }

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

  417.                 // this was the last point we needed, we can compute the derivatives
  418.                 stepStart = t[0];
  419.                 stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);

  421.                 // first scaled derivative
  422.                 scaled = yDot[0].clone();
  423.                 for (int j = 0; j < scaled.length; ++j) {
  424.                     scaled[j] *= stepSize;
  425.                 }

  427.                 // higher order derivatives
  428.                 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);

  430.                 // stop the integrator now that all needed steps have been handled
  431.                 throw new InitializationCompletedMarkerException();
  432.             }
  433.         }

  435.         /** {@inheritDoc} */
  436.         @Override
  437.         public void init(double t0, double[] y0, double time) {
  438.             // nothing to do
  439.         }
  440.     }

  442.     /** Marker exception used ONLY to stop the starter integrator after first step. */
  443.     private static final class InitializationCompletedMarkerException
  444.         extends RuntimeException {

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

  449.         /** Simple constructor. */
  450.         InitializationCompletedMarkerException() {
  451.             super((Throwable) null);
  452.         }
  453.     }
  454. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT