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.  *      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.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.MathIllegalStateException;
  24. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  25. import org.apache.commons.math4.legacy.exception.NoBracketingException;
  26. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  27. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  28. import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix;
  29. import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
  30. import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853FieldIntegrator;
  31. import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
  32. import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
  33. import org.apache.commons.math4.core.jdkmath.JdkMath;
  34. import org.apache.commons.math4.legacy.core.MathArrays;

  36. /**
  37.  * This class is the base class for multistep integrators for Ordinary
  38.  * Differential Equations.
  39.  * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
  40.  * <div style="white-space: pre"><code>
  41.  * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
  42.  * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
  43.  * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
  44.  * ...
  45.  * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
  46.  * </code></div>
  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.  * <div style="white-space: pre"><code>
  51.  * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
  52.  * </code></div>
  53.  * (we omit the k index in the notation for clarity)
  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.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthFieldIntegrator
  63.  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonFieldIntegrator
  64.  * @param <T> the type of the field elements
  65.  * @since 3.6
  66.  */
  67. public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
  68.     extends AdaptiveStepsizeFieldIntegrator<T> {

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

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

  78.     /** Starter integrator. */
  79.     private FirstOrderFieldIntegrator<T> starter;

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

  84.     /** Stepsize control exponent. */
  85.     private double exp;

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

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

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

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

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

  126.         if (nSteps < 2) {
  127.             throw new NumberIsTooSmallException(
  128.                   LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
  129.                   nSteps, 2, true);
  130.         }

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

  137.         exp = -1.0 / order;

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

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

  176.         exp = -1.0 / order;

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

  184.     /**
  185.      * Get the starter integrator.
  186.      * @return starter integrator
  187.      */
  188.     public FirstOrderFieldIntegrator<T> getStarterIntegrator() {
  189.         return starter;
  190.     }

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

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

  225.         // make sure NO user event nor user step handler is triggered,
  226.         // this is the task of the top level integrator, not the task
  227.         // of the starter integrator
  228.         starter.clearEventHandlers();
  229.         starter.clearStepHandlers();

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

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

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

  239.             // we should not reach this step
  240.             throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
  241.         } catch (InitializationCompletedMarkerException icme) { // NOPMD
  242.             // this is the expected nominal interruption of the start integrator

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

  248.         // remove the specific step handler
  249.         starter.clearStepHandlers();
  250.     }

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

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

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

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

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

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

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

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

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

  320.         final T ratio = newStepSize.divide(getStepSize());
  321.         for (int i = 0; i < scaled.length; ++i) {
  322.             scaled[i] = scaled[i].multiply(ratio);
  323.         }

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

  335.         setStepSize(newStepSize);
  336.     }


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

  349.     /** Specialized step handler storing the first step.
  350.      */
  351.     private final class FieldNordsieckInitializer implements FieldStepHandler<T> {

  353.         /** Equation mapper. */
  354.         private final FieldEquationsMapper<T> mapper;

  356.         /** Steps counter. */
  357.         private int count;

  359.         /** Saved start. */
  360.         private FieldODEStateAndDerivative<T> savedStart;

  362.         /** First steps times. */
  363.         private final T[] t;

  365.         /** First steps states. */
  366.         private final T[][] y;

  368.         /** First steps derivatives. */
  369.         private final T[][] yDot;

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

  383.         /** {@inheritDoc} */
  384.         @Override
  385.         public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
  386.             throws MaxCountExceededException {


  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]    = mapper.mapState(prev);
  395.                 yDot[count] = mapper.mapDerivative(prev);
  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]    = mapper.mapState(curr);
  403.             yDot[count] = mapper.mapDerivative(curr);

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

  407.                 // this was the last point we needed, we can compute the derivatives
  408.                 setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));

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

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

  419.                 // stop the integrator now that all needed steps have been handled
  420.                 setStepStart(savedStart);
  421.                 throw new InitializationCompletedMarkerException();
  422.             }
  423.         }

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

  432.     /** Marker exception used ONLY to stop the starter integrator after first step. */
  433.     private static final class InitializationCompletedMarkerException
  434.         extends RuntimeException {

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

  439.         /** Simple constructor. */
  440.         InitializationCompletedMarkerException() {
  441.             super((Throwable) null);
  442.         }
  443.     }
  444. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT