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:
|
| 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. |
| 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 final 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. |
| Access | Constructor and Description |
|---|---|
| public | GraggBulirschStoerIntegrator(final double
minimal step (sign is irrelevant, regardless of
integration direction, forward or backward), the last step can
be smaller than this minStep, final double maximal step (sign is irrelevant, regardless of
integration direction, forward or backward), the last step can
be smaller than this maxStep, final double allowed absolute error scalAbsoluteTolerance, final double allowed relative error scalRelativeTolerance)Simple constructor. |
| public | GraggBulirschStoerIntegrator(final double
minimal step (must be positive even for backward
integration), the last step can be smaller than this minStep, final double maximal step (must be positive even for backward
integration) maxStep, final double[] allowed absolute error vecAbsoluteTolerance, final double[] allowed relative error vecRelativeTolerance)Simple constructor. |
| Modifier and Type | Method and Description |
|---|---|
| public void | addEventHandler(final EventHandler
event handler function, final double maximal time interval between switching
function checks (this interval prevents missing sign changes in
case the integration steps becomes very large) maxCheckInterval, final double convergence threshold in the event time search convergence, final int upper limit of the iteration count in
the event time search maxIterationCount, final UnivariateSolver The root-finding algorithm to use to detect the state
events. solver)Overrides org. Implements org. Add an event handler to the integrator. |
| public void | addStepHandler(final StepHandler
handler for the accepted steps handler)Overrides org. Implements org. Add a step handler to this integrator. |
| private void | extrapolate(final int
offset to use in the coefficients table offset, final int index of the last updated point k, final double[][] working diagonal of the Aitken-Neville's
triangle, without the last element diag, final double[] last element last)Extrapolate a vector. |
| private void | |
| public void | integrate(final ExpandableStatefulODE
complete set of differential equations to integrate equations, final double target time for the integration
(can be set to a value smaller than t)t0 for backward integration)Implements abstract org. Integrate a set of differential equations up to the given time. |
| private void | rescale(final double[]
first state vector to use for scaling y1, final double[] second state vector to use for scaling y2, final double[] scaling array to update (can be shorter than state) scale)Update scaling array. |
| public void | setControlFactors(final double
first stepsize control factor (the factor is
reset to default if lower than 0.0001 or greater than 0.9999) control1, final double second stepsize control factor (the factor
is reset to default if lower than 0.0001 or greater than 0.9999) control2, final double third stepsize control factor (the factor is
reset to default if lower than 0.0001 or greater than 0.9999) control3, final double fourth stepsize control factor (the factor
is reset to default if lower than 1.0001 or greater than 999.9) control4)Set the step size control factors. |
| public void | setInterpolationControl(final boolean
if true, interpolation error is used
for stepsize control useInterpolationErrorForControl, final int interpolation order control parameter (the parameter
is reset to default if <= 0 or >= 7) mudifControlParameter)Set the interpolation order control parameter. |
| public void | setOrderControl(final int
maximal order in the extrapolation table (the
maximal order is reset to default if order <= 6 or odd) maximalOrder, final double first order control factor (the factor is
reset to default if lower than 0.0001 or greater than 0.9999) control1, final double second order control factor (the factor
is reset to default if lower than 0.0001 or greater than 0.9999) control2)Set the order control parameters. |
| public void | setStabilityCheck(final boolean
if true, stability check will be performed,
if false, the check will be skipped performStabilityCheck, final int maximal number of iterations for which checks are
performed (the number of iterations is reset to default if negative
or null) maxNumIter, final int maximal number of checks for each iteration
(the number of checks is reset to default if negative or null) maxNumChecks, final double stepsize reduction factor in case of
failure (the factor is reset to default if lower than 0.0001 or
greater than 0.9999) stepsizeReductionFactor)Set the stability check controls. |
| private boolean | Returns: true if computation was done properly, false if stability check failed before end of computationinitial time t0, final double[] initial value of the state vector at t0 y0, final double global step step, final int iteration number (from 0 to sequence.length - 1) k, final double[] scaling array (can be shorter than state) scale, final double[][] placeholder where to put the state vector derivatives at each substep
(element 0 already contains initial derivative) f, final double[] placeholder where to put the state vector at the middle of the step yMiddle, final double[] placeholder where to put the state vector at the end yEnd, final double[] placeholder for one state vector yTmp)Perform integration over one step using substeps of a modified midpoint method. |
| coeff | back to summary |
|---|---|
| private double[][] coeff extrapolation coefficients. | |
| costPerStep | back to summary |
|---|---|
| private int[] costPerStep overall cost of applying step reduction up to iteration k+1, in number of calls. | |
| costPerTimeUnit | back to summary |
|---|---|
| private double[] costPerTimeUnit cost per unit step. | |
| maxChecks | back to summary |
|---|---|
| private int maxChecks maximal number of checks for each iteration. | |
| maxIter | back to summary |
|---|---|
| private int maxIter maximal number of iterations for which checks are performed. | |
| maxOrder | back to summary |
|---|---|
| private int maxOrder maximal order. | |
| METHOD_NAME | back to summary |
|---|---|
| private static final String METHOD_NAME Integrator method name. | |
| mudif | back to summary |
|---|---|
| private int mudif interpolation order control parameter. | |
| optimalStep | back to summary |
|---|---|
| private double[] optimalStep optimal steps for each order. | |
| orderControl1 | back to summary |
|---|---|
| private double orderControl1 first order control factor. | |
| orderControl2 | back to summary |
|---|---|
| private double orderControl2 second order control factor. | |
| performTest | back to summary |
|---|---|
| private boolean performTest stability check enabling parameter. | |
| sequence | back to summary |
|---|---|
| private int[] sequence step size sequence. | |
| stabilityReduction | back to summary |
|---|---|
| private double stabilityReduction stepsize reduction factor in case of stability check failure. | |
| stepControl1 | back to summary |
|---|---|
| private double stepControl1 first stepsize control factor. | |
| stepControl2 | back to summary |
|---|---|
| private double stepControl2 second stepsize control factor. | |
| stepControl3 | back to summary |
|---|---|
| private double stepControl3 third stepsize control factor. | |
| stepControl4 | back to summary |
|---|---|
| private double stepControl4 fourth stepsize control factor. | |
| useInterpolationError | back to summary |
|---|---|
| private boolean useInterpolationError use interpolation error in stepsize control. | |
| GraggBulirschStoerIntegrator | back to summary |
|---|---|
| public GraggBulirschStoerIntegrator(final double minStep, final double maxStep, final double scalAbsoluteTolerance, final 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.
| |
| GraggBulirschStoerIntegrator | back to summary |
|---|---|
| public GraggBulirschStoerIntegrator(final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final 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.
| |
| addEventHandler | back to summary |
|---|---|
| public void addEventHandler(final EventHandler function, final double maxCheckInterval, final double convergence, final int maxIterationCount, final UnivariateSolver solver) Overrides org. Implements org. Doc from org. Add an event handler to the integrator.
| |
| addStepHandler | back to summary |
|---|---|
| public void addStepHandler(final StepHandler handler) Overrides org. Implements org. Doc from org. Add a step handler to this integrator. The handler will be called by the integrator for each accepted step.
| |
| extrapolate | back to summary |
|---|---|
| private void extrapolate(final int offset, final int k, final double[][] diag, final double[] last) Extrapolate a vector.
| |
| initializeArrays | back to summary |
|---|---|
| private void initializeArrays() Initialize the integrator internal arrays. | |
| integrate | back to summary |
|---|---|
| public void integrate(final ExpandableStatefulODE equations, final double t) throws NumberIsTooSmallException, DimensionMismatchException, MaxCountExceededException, NoBracketingException Implements abstract org. Doc from org. 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 (
| |
| rescale | back to summary |
|---|---|
| private void rescale(final double[] y1, final double[] y2, final double[] scale) Update scaling array.
| |
| setControlFactors | back to summary |
|---|---|
| public void setControlFactors(final double control1, final double control2, final double control3, final 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.
| |
| setInterpolationControl | back to summary |
|---|---|
| public void setInterpolationControl(final boolean useInterpolationErrorForControl, final 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.
| |
| setOrderControl | back to summary |
|---|---|
| public void setOrderControl(final int maximalOrder, final double control1, final 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.
| |
| setStabilityCheck | back to summary |
|---|---|
| public void setStabilityCheck(final boolean performStabilityCheck, final int maxNumIter, final int maxNumChecks, final 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.
| |
| tryStep | back to summary |
|---|---|
| private boolean tryStep(final double t0, final double[] y0, final double step, final int k, final double[] scale, final double[][] f, final double[] yMiddle, final double[] yEnd, final double[] yTmp) throws MaxCountExceededException, DimensionMismatchException Perform integration over one step using substeps of a modified midpoint method.
| |