org.jscience.mathematics.analysis.ode
Class GraggBulirschStoerIntegrator

java.lang.Object
  extended by org.jscience.mathematics.analysis.ode.AdaptiveStepsizeIntegrator
      extended by org.jscience.mathematics.analysis.ode.GraggBulirschStoerIntegrator
All Implemented Interfaces:
FirstOrderIntegrator, Named

public class GraggBulirschStoerIntegrator
extends AdaptiveStepsizeIntegrator
implements Named

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 Runge-Kutta-Fehlberg 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.


Field Summary
 
Fields inherited from class org.jscience.mathematics.analysis.ode.AdaptiveStepsizeIntegrator
handler, scalAbsoluteTolerance, scalRelativeTolerance, switchesHandler, vecAbsoluteTolerance, vecRelativeTolerance
 
Constructor Summary
GraggBulirschStoerIntegrator(double minStep, double maxStep, double[] vecAbsoluteTolerance, double[] vecRelativeTolerance)
          Simple constructor.
GraggBulirschStoerIntegrator(double minStep, double maxStep, double scalAbsoluteTolerance, double scalRelativeTolerance)
          Simple constructor.
 
Method Summary
 void addSwitchingFunction(SwitchingFunction function, double maxCheckInterval, double convergence)
          Add a switching function to the integrator.
 java.lang.String getName()
          Get the name of the method.
 void integrate(FirstOrderDifferentialEquations equations, double t0, double[] y0, double t, double[] y)
          Integrate the differential equations up to the given time
 void setInterpolationControl(boolean useInterpolationError, int mudif)
          Set the interpolation order control parameter.
 void setOrderControl(int maxOrder, double orderControl1, double orderControl2)
          Set the order control parameters.
 void setStabilityCheck(boolean performTest, int maxIter, int maxChecks, double stabilityReduction)
          Set the stability check controls.
 void setStepHandler(StepHandler handler)
          Set the step handler for this integrator.
 void setStepsizeControl(double stepControl1, double stepControl2, double stepControl3, double stepControl4)
          Set the step size control factors.
 
Methods inherited from class org.jscience.mathematics.analysis.ode.AdaptiveStepsizeIntegrator
filterStep, getMaxStep, getMinStep, getStepHandler, initializeStep, setInitialStepSize
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

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 (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)
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 performTest,
                              int maxIter,
                              int maxChecks,
                              double stabilityReduction)
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:
performTest - if true, stability check will be performed, if false, the check will be skipped
maxIter - maximal number of iterations for which checks are performed (the number of iterations is reset to default if negative or null)
maxChecks - maximal number of checks for each iteration (the number of checks is reset to default if negative or null)
stabilityReduction - stepsize reduction factor in case of failure (the factor is reset to default if lower than 0.0001 or greater than 0.9999)

setStepsizeControl

public void setStepsizeControl(double stepControl1,
                               double stepControl2,
                               double stepControl3,
                               double stepControl4)
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:
stepControl1 - first stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
stepControl2 - second stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
stepControl3 - third stepsize control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
stepControl4 - 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 maxOrder,
                            double orderControl1,
                            double orderControl2)
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:
maxOrder - maximal order in the extrapolation table (the maximal order is reset to default if order <= 6 or odd)
orderControl1 - first order control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)
orderControl2 - second order control factor (the factor is reset to default if lower than 0.0001 or greater than 0.9999)

setStepHandler

public void setStepHandler(StepHandler handler)
Set the step handler for this integrator. The handler will be called by the integrator for each accepted step.

Specified by:
setStepHandler in interface FirstOrderIntegrator
Overrides:
setStepHandler in class AdaptiveStepsizeIntegrator
Parameters:
handler - handler for the accepted steps

addSwitchingFunction

public void addSwitchingFunction(SwitchingFunction function,
                                 double maxCheckInterval,
                                 double convergence)
Add a switching function to the integrator.

Specified by:
addSwitchingFunction in interface FirstOrderIntegrator
Overrides:
addSwitchingFunction in class AdaptiveStepsizeIntegrator
Parameters:
function - switching function
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

setInterpolationControl

public void setInterpolationControl(boolean useInterpolationError,
                                    int mudif)
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:
useInterpolationError - if true, interpolation error is used for stepsize control
mudif - interpolation order control parameter (the parameter is reset to default if <= 0 or >= 7)

getName

public java.lang.String getName()
Get the name of the method.

Specified by:
getName in interface FirstOrderIntegrator
Specified by:
getName in interface Named
Returns:
name of the method

integrate

public void integrate(FirstOrderDifferentialEquations equations,
                      double t0,
                      double[] y0,
                      double t,
                      double[] y)
               throws DerivativeException,
                      IntegrationException
Description copied from interface: FirstOrderIntegrator
Integrate the differential equations up to the given time

Specified by:
integrate in interface FirstOrderIntegrator
Specified by:
integrate in class AdaptiveStepsizeIntegrator
Parameters:
equations - differential equations to integrate
t0 - initial time
y0 - initial value of the state vector at t0
t - target time for the integration (can be set to a value smaller thant t0 for backward integration)
y - placeholder where to put the state vector at each successful step (and hence at the end of integration), can be the same object as y0
Throws:
DerivativeException - this exception is propagated to the caller if the underlying user function triggers one
IntegrationException - if the integrator cannot perform integration