org.apache.commons.math3.ode.nonstiff

## Class GraggBulirschStoerIntegrator

- java.lang.Object
- - org.apache.commons.math3.ode.AbstractIntegrator
  - - org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator
    - - org.apache.commons.math3.ode.nonstiff.GraggBulirschStoerIntegrator

- All Implemented Interfaces:
    FirstOrderIntegrator, ODEIntegrator

      public class GraggBulirschStoerIntegrator
      extends AdaptiveStepsizeIntegrator

  This class implements a Gragg-Bulirsch-Stoer integrator for Ordinary Differential Equations.

  The Gragg-Bulirsch-Stoer algorithm is one of the most efficient ones currently available for smooth problems. It uses Richardson extrapolation to estimate what would be the solution if the step size could be decreased down to zero.

  This method changes both the step size and the order during integration, in order to minimize computation cost. It is particularly well suited when a very high precision is needed. The limit where this method becomes more efficient than high-order embedded Runge-Kutta methods like `Dormand-Prince 8(5,3)` depends on the problem. Results given in the Hairer, Norsett and Wanner book show for example that this limit occurs for accuracy around 1e-6 when integrating Saltzam-Lorenz equations (the authors note this problem is extremely sensitive to the errors in the first integration steps), and around 1e-11 for a two dimensional celestial mechanics problems with seven bodies (pleiades problem, involving quasi-collisions for which automatic step size control is essential).

  This implementation is basically a reimplementation in Java of the odex fortran code by E. Hairer and G. Wanner. The redistribution policy for this code is available here, for convenience, it is reproduced below.

  Copyright (c) 2004, Ernst Hairer
  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

  Since:
    1.2
  Version:
    $Id: GraggBulirschStoerIntegrator.java 1416643 2012-12-03 19:37:14Z tn $

- - ### Field Summary

    Fields
    Modifier and Type | Field and Description
    --- | ---
    `private double[][]` | `coeff` extrapolation coefficients.
    `private int[]` | `costPerStep` overall cost of applying step reduction up to iteration k+1, in number of calls.
    `private double[]` | `costPerTimeUnit` cost per unit step.
    `private int` | `maxChecks` maximal number of checks for each iteration.
    `private int` | `maxIter` maximal number of iterations for which checks are performed.
    `private int` | `maxOrder` maximal order.
    `private static String` | `METHOD_NAME` Integrator method name.
    `private int` | `mudif` interpolation order control parameter.
    `private double[]` | `optimalStep` optimal steps for each order.
    `private double` | `orderControl1` first order control factor.
    `private double` | `orderControl2` second order control factor.
    `private boolean` | `performTest` stability check enabling parameter.
    `private int[]` | `sequence` step size sequence.
    `private double` | `stabilityReduction` stepsize reduction factor in case of stability check failure.
    `private double` | `stepControl1` first stepsize control factor.
    `private double` | `stepControl2` second stepsize control factor.
    `private double` | `stepControl3` third stepsize control factor.
    `private double` | `stepControl4` fourth stepsize control factor.
    `private boolean` | `useInterpolationError` use interpolation error in stepsize control.

    - ### Fields inherited from class org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator

      `mainSetDimension, scalAbsoluteTolerance, scalRelativeTolerance, vecAbsoluteTolerance, vecRelativeTolerance`

    - ### Fields inherited from class org.apache.commons.math3.ode.AbstractIntegrator

      `isLastStep, resetOccurred, stepHandlers, stepSize, stepStart`

  - ### Constructor Summary

    Constructors
    Constructor and Description
    `GraggBulirschStoerIntegrator(double minStep, double maxStep, double[] vecAbsoluteTolerance, double[] vecRelativeTolerance)` Simple constructor.
    `GraggBulirschStoerIntegrator(double minStep, double maxStep, double scalAbsoluteTolerance, double scalRelativeTolerance)` Simple constructor.

  - ### Method Summary

    Methods
    Modifier and Type | Method and Description
    --- | ---
    `void` | `addEventHandler(EventHandler function, double maxCheckInterval, double convergence, int maxIterationCount, UnivariateSolver solver)` Add an event handler to the integrator.
    `void` | `addStepHandler(StepHandler handler)` Add a step handler to this integrator.
    `private void` | `extrapolate(int offset, int k, double[][] diag, double[] last)` Extrapolate a vector.
    `private void` | `initializeArrays()` Initialize the integrator internal arrays.
    `void` | `integrate(ExpandableStatefulODE equations, double t)` Integrate a set of differential equations up to the given time.
    `private void` | `rescale(double[] y1, double[] y2, double[] scale)` Update scaling array.
    `void` | `setControlFactors(double control1, double control2, double control3, double control4)` Set the step size control factors.
    `void` | `setInterpolationControl(boolean useInterpolationErrorForControl, int mudifControlParameter)` Set the interpolation order control parameter.
    `void` | `setOrderControl(int maximalOrder, double control1, double control2)` Set the order control parameters.
    `void` | `setStabilityCheck(boolean performStabilityCheck, int maxNumIter, int maxNumChecks, double stepsizeReductionFactor)` Set the stability check controls.
    `private boolean` | `tryStep(double t0, double[] y0, double step, int k, double[] scale, double[][] f, double[] yMiddle, double[] yEnd, double[] yTmp)` Perform integration over one step using substeps of a modified midpoint method.

    - ### Methods inherited from class org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator

      `filterStep, getCurrentStepStart, getMaxStep, getMinStep, initializeStep, resetInternalState, sanityChecks, setInitialStepSize, setStepSizeControl, setStepSizeControl`

    - ### Methods inherited from class org.apache.commons.math3.ode.AbstractIntegrator

      `acceptStep, addEventHandler, clearEventHandlers, clearStepHandlers, computeDerivatives, getCurrentSignedStepsize, getEvaluations, getEventHandlers, getMaxEvaluations, getName, getStepHandlers, initIntegration, integrate, setEquations, setMaxEvaluations, setStateInitialized`

    - ### Methods inherited from class java.lang.Object

      `clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait`

- - ### Field Detail

    - #### METHOD_NAME

          private static final String METHOD_NAME

      Integrator method name.
      See Also:
        Constant Field Values

    - #### maxOrder

          private int maxOrder

      maximal order.

    - #### sequence

          private int[] sequence

      step size sequence.

    - #### costPerStep

          private int[] costPerStep

      overall cost of applying step reduction up to iteration k+1, in number of calls.

    - #### costPerTimeUnit

          private double[] costPerTimeUnit

      cost per unit step.

    - #### optimalStep

          private double[] optimalStep

      optimal steps for each order.

    - #### coeff

          private double[][] coeff

      extrapolation coefficients.

    - #### performTest

          private boolean performTest

      stability check enabling parameter.

    - #### maxChecks

          private int maxChecks

      maximal number of checks for each iteration.

    - #### maxIter

          private int maxIter

      maximal number of iterations for which checks are performed.

    - #### stabilityReduction

          private double stabilityReduction

      stepsize reduction factor in case of stability check failure.

    - #### stepControl1

          private double stepControl1

      first stepsize control factor.

    - #### stepControl2

          private double stepControl2

      second stepsize control factor.

    - #### stepControl3

          private double stepControl3

      third stepsize control factor.

    - #### stepControl4

          private double stepControl4

      fourth stepsize control factor.

    - #### orderControl1

          private double orderControl1

      first order control factor.

    - #### orderControl2

          private double orderControl2

      second order control factor.

    - #### useInterpolationError

          private boolean useInterpolationError

      use interpolation error in stepsize control.

    - #### mudif

          private int mudif

      interpolation order control parameter.

  - ### Constructor Detail

    - #### GraggBulirschStoerIntegrator

          public GraggBulirschStoerIntegrator(double minStep,
                                      double maxStep,
                                      double scalAbsoluteTolerance,
                                      double scalRelativeTolerance)

      Simple constructor. Build a Gragg-Bulirsch-Stoer integrator with the given step bounds. All tuning parameters are set to their default values. The default step handler does nothing.
      Parameters:
        `minStep` - minimal step (sign is irrelevant, regardless of integration direction, forward or backward), the last step can be smaller than this
        `maxStep` - maximal step (sign is irrelevant, regardless of integration direction, forward or backward), the last step can be smaller than this
        `scalAbsoluteTolerance` - allowed absolute error
        `scalRelativeTolerance` - allowed relative error

    - #### GraggBulirschStoerIntegrator

          public GraggBulirschStoerIntegrator(double minStep,
                                      double maxStep,
                                      double[] vecAbsoluteTolerance,
                                      double[] vecRelativeTolerance)

      Simple constructor. Build a Gragg-Bulirsch-Stoer integrator with the given step bounds. All tuning parameters are set to their default values. The default step handler does nothing.
      Parameters:
        `minStep` - minimal step (must be positive even for backward integration), the last step can be smaller than this
        `maxStep` - maximal step (must be positive even for backward integration)
        `vecAbsoluteTolerance` - allowed absolute error
        `vecRelativeTolerance` - allowed relative error

  - ### Method Detail

    - #### setStabilityCheck

          public void setStabilityCheck(boolean performStabilityCheck,
                               int maxNumIter,
                               int maxNumChecks,
                               double stepsizeReductionFactor)

      Set the stability check controls.

      The stability check is performed on the first few iterations of the extrapolation scheme. If this test fails, the step is rejected and the stepsize is reduced.

      By default, the test is performed, at most during two iterations at each step, and at most once for each of these iterations. The default stepsize reduction factor is 0.5.

      Parameters:
        `performStabilityCheck` - if true, stability check will be performed, if false, the check will be skipped
        `maxNumIter` - maximal number of iterations for which checks are performed (the number of iterations is reset to default if negative or null)
        `maxNumChecks` - maximal number of checks for each iteration (the number of checks is reset to default if negative or null)
        `stepsizeReductionFactor` - stepsize reduction factor in case of failure (the factor is reset to default if lower than 0.0001 or greater than 0.9999)

    - #### setControlFactors

          public void setControlFactors(double control1,
                               double control2,
                               double control3,
                               double control4)

      Set the step size control factors.

      The new step size hNew is computed from the old one h by:

           hNew = h * stepControl2 / (err/stepControl1)^(1/(2k+1))

      where err is the scaled error and k the iteration number of the extrapolation scheme (counting from 0). The default values are 0.65 for stepControl1 and 0.94 for stepControl2.

      The step size is subject to the restriction:

           stepControl3^(1/(2k+1))/stepControl4 <= hNew/h <= 1/stepControl3^(1/(2k+1))

      The default values are 0.02 for stepControl3 and 4.0 for stepControl4.
      Parameters:
        `control1` - first stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
        `control2` - second stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
        `control3` - third stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
        `control4` - fourth stepsize control factor (the factor is reset to default if lower than 1.0001 or greater than 999.9)

    - #### setOrderControl

          public void setOrderControl(int maximalOrder,
                             double control1,
                             double control2)

      Set the order control parameters.

      The Gragg-Bulirsch-Stoer method changes both the step size and the order during integration, in order to minimize computation cost. Each extrapolation step increases the order by 2, so the maximal order that will be used is always even, it is twice the maximal number of columns in the extrapolation table.

           order is decreased if w(k-1) <= w(k)   * orderControl1
           order is increased if w(k)   <= w(k-1) * orderControl2

      where w is the table of work per unit step for each order (number of function calls divided by the step length), and k is the current order.

      The default maximal order after construction is 18 (i.e. the maximal number of columns is 9). The default values are 0.8 for orderControl1 and 0.9 for orderControl2.

      Parameters:
        `maximalOrder` - maximal order in the extrapolation table (the maximal order is reset to default if order <= 6 or odd)
        `control1` - first order control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
        `control2` - second order control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)

    - #### addStepHandler

          public void addStepHandler(StepHandler handler)

      Add a step handler to this integrator.

      The handler will be called by the integrator for each accepted step.

      Specified by:
        `addStepHandler` in interface `ODEIntegrator`
      Overrides:
        `addStepHandler` in class `AbstractIntegrator`
      Parameters:
        `handler` - handler for the accepted steps
      See Also:
        `ODEIntegrator.getStepHandlers()`, `ODEIntegrator.clearStepHandlers()`

    - #### addEventHandler

          public void addEventHandler(EventHandler function,
                             double maxCheckInterval,
                             double convergence,
                             int maxIterationCount,
                             UnivariateSolver solver)

      Add an event handler to the integrator.
      Specified by:
        `addEventHandler` in interface `ODEIntegrator`
      Overrides:
        `addEventHandler` in class `AbstractIntegrator`
      Parameters:
        `function` - event handler
        `maxCheckInterval` - maximal time interval between switching function checks (this interval prevents missing sign changes in case the integration steps becomes very large)
        `convergence` - convergence threshold in the event time search
        `maxIterationCount` - upper limit of the iteration count in the event time search
        `solver` - The root-finding algorithm to use to detect the state events.
      See Also:
        `ODEIntegrator.getEventHandlers()`, `ODEIntegrator.clearEventHandlers()`

    - #### initializeArrays

          private void initializeArrays()

      Initialize the integrator internal arrays.

    - #### setInterpolationControl

          public void setInterpolationControl(boolean useInterpolationErrorForControl,
                                     int mudifControlParameter)

      Set the interpolation order control parameter. The interpolation order for dense output is 2k - mudif + 1. The default value for mudif is 4 and the interpolation error is used in stepsize control by default.
      Parameters:
        `useInterpolationErrorForControl` - if true, interpolation error is used for stepsize control
        `mudifControlParameter` - interpolation order control parameter (the parameter is reset to default if <= 0 or >= 7)

    - #### rescale

          private void rescale(double[] y1,
                     double[] y2,
                     double[] scale)

      Update scaling array.
      Parameters:
        `y1` - first state vector to use for scaling
        `y2` - second state vector to use for scaling
        `scale` - scaling array to update (can be shorter than state)

    - #### tryStep

          private boolean tryStep(double t0,
                        double[] y0,
                        double step,
                        int k,
                        double[] scale,
                        double[][] f,
                        double[] yMiddle,
                        double[] yEnd,
                        double[] yTmp)
                           throws MaxCountExceededException,
                                  DimensionMismatchException

      Perform integration over one step using substeps of a modified midpoint method.
      Parameters:
        `t0` - initial time
        `y0` - initial value of the state vector at t0
        `step` - global step
        `k` - iteration number (from 0 to sequence.length - 1)
        `scale` - scaling array (can be shorter than state)
        `f` - placeholder where to put the state vector derivatives at each substep (element 0 already contains initial derivative)
        `yMiddle` - placeholder where to put the state vector at the middle of the step
        `yEnd` - placeholder where to put the state vector at the end
        `yTmp` - placeholder for one state vector
      Returns:
        true if computation was done properly, false if stability check failed before end of computation
      Throws:
        `MaxCountExceededException` - if the number of functions evaluations is exceeded
        `DimensionMismatchException` - if arrays dimensions do not match equations settings

    - #### extrapolate

          private void extrapolate(int offset,
                         int k,
                         double[][] diag,
                         double[] last)

      Extrapolate a vector.
      Parameters:
        `offset` - offset to use in the coefficients table
        `k` - index of the last updated point
        `diag` - working diagonal of the Aitken-Neville's triangle, without the last element
        `last` - last element

    - #### integrate

          public void integrate(ExpandableStatefulODE equations,
                       double t)
                         throws NumberIsTooSmallException,
                                DimensionMismatchException,
                                MaxCountExceededException,
                                NoBracketingException

      Integrate a set of differential equations up to the given time.

      This method solves an Initial Value Problem (IVP).

      The set of differential equations is composed of a main set, which can be extended by some sets of secondary equations. The set of equations must be already set up with initial time and partial states. At integration completion, the final time and partial states will be available in the same object.

      Since this method stores some internal state variables made available in its public interface during integration (`AbstractIntegrator.getCurrentSignedStepsize()`), it is not thread-safe.

      Specified by:
        `integrate` in class `AdaptiveStepsizeIntegrator`
      Parameters:
        `equations` - complete set of differential equations to integrate
        `t` - target time for the integration (can be set to a value smaller than `t0` for backward integration)
      Throws:
        `NumberIsTooSmallException` - if integration step is too small
        `DimensionMismatchException` - if the dimension of the complete state does not match the complete equations sets dimension
        `MaxCountExceededException` - if the number of functions evaluations is exceeded
        `NoBracketingException` - if the location of an event cannot be bracketed

Copyright (c) 2003-2013 Apache Software Foundation