Go to the first, previous, next, last section, table of contents.
This chapter describes functions for solving ordinary differential
equation (ODE) initial value problems. The library provides a variety
of lowlevel methods, such as RungeKutta and BulirschStoer routines,
and higherlevel components for adaptive stepsize control. The
components can be combined by the user to achieve the desired
solution, with full access to any intermediate steps. A driver object
can be used as a high level wrapper for easy use of low level
functions.
These functions are declared in the header file `gsl_odeiv2.h'.
This is a new interface in version 1.15 and uses the prefix
gsl_odeiv2
for all functions. It is recommended over the
previous gsl_odeiv
implementation defined in `gsl_odeiv.h'
The old interface has been retained under the original name for
backwards compatibility.
The routines solve the general ndimensional firstorder system,
dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t))
for i = 1, \dots, n. The stepping functions rely on the vector
of derivatives f_i and the Jacobian matrix,
J_{ij} = df_i(t,y(t)) / dy_j.
A system of equations is defined using the gsl_odeiv2_system
datatype.
 Data Type: gsl_odeiv2_system

This data type defines a general ODE system with arbitrary parameters.
int (* function) (double t, const double y[], double dydt[], void * params)

This function should store the vector elements
f_i(t,y,params) in the array dydt,
for arguments (t,y) and parameters params.
The function should return
GSL_SUCCESS
if the calculation was
completed successfully. Any other return value indicates an error. A
special return value GSL_EBADFUNC
causes gsl_odeiv2
routines to immediately stop and return. The user must call an
appropriate reset function (e.g. gsl_odeiv2_driver_reset
or
gsl_odeiv2_step_reset
) before continuing. Use return values
distinct from standard GSL error codes to distinguish your function as
the source of the error.
int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);

This function should store the vector of derivative elements
in the array dfdt and the Jacobian matrix
J_{ij} in the array dfdy, regarded as a rowordered
matrix
J(i,j) = dfdy[i * dimension + j]
where dimension
is the dimension of the system.
Not all of the stepper algorithms of gsl_odeiv2
make use of the
Jacobian matrix, so it may not be necessary to provide this function
(the jacobian
element of the struct can be replaced by a null
pointer for those algorithms).
The function should return GSL_SUCCESS
if the calculation was
completed successfully. Any other return value indicates an error. A
special return value GSL_EBADFUNC
causes gsl_odeiv2
routines to immediately stop and return. The user must call an
appropriate reset function (e.g. gsl_odeiv2_driver_reset
or
gsl_odeiv2_step_reset
) before continuing. Use return values
distinct from standard GSL error codes to distinguish your function as
the source of the error.
size_t dimension;

This is the dimension of the system of equations.
void * params

This is a pointer to the arbitrary parameters of the system.
The lowest level components are the stepping functions which
advance a solution from time t to t+h for a fixed
stepsize h and estimate the resulting local error.
 Function: gsl_odeiv2_step * gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T, size_t dim)

This function returns a pointer to a newly allocated instance of a
stepping function of type T for a system of dim
dimensions. Please note that if you use a stepper method that
requires access to a driver object, it is advisable to use a driver
allocation method, which automatically allocates a stepper, too.
 Function: int gsl_odeiv2_step_reset (gsl_odeiv2_step * s)

This function resets the stepping function s. It should be used
whenever the next use of s will not be a continuation of a
previous step.
 Function: void gsl_odeiv2_step_free (gsl_odeiv2_step * s)

This function frees all the memory associated with the stepping function
s.
 Function: const char * gsl_odeiv2_step_name (const gsl_odeiv2_step * s)

This function returns a pointer to the name of the stepping function.
For example,
printf ("step method is '%s'\n",
gsl_odeiv2_step_name (s));
would print something like step method is 'rkf45'
.
 Function: unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * s)

This function returns the order of the stepping function on the previous
step. The order can vary if the stepping function itself is adaptive.
 Function: int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s, const gsl_odeiv2_driver * d)

This function sets a pointer of the driver object d for stepper
s, to allow the stepper to access control (and evolve) object
through the driver object. This is a requirement for some steppers, to
get the desired error level for internal iteration of
stepper. Allocation of a driver object calls this function
automatically.
 Function: int gsl_odeiv2_step_apply (gsl_odeiv2_step * s, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv2_system * sys)

This function applies the stepping function s to the system of
equations defined by sys, using the stepsize h to advance
the system from time t and state y to time t+h.
The new state of the system is stored in y on output, with an
estimate of the absolute error in each component stored in yerr.
If the argument dydt_in is not null it should point an array
containing the derivatives for the system at time t on input. This
is optional as the derivatives will be computed internally if they are
not provided, but allows the reuse of existing derivative information.
On output the new derivatives of the system at time t+h will
be stored in dydt_out if it is not null.
The stepping function returns GSL_FAILURE
if it is unable to
compute the requested step. Also, if the usersupplied functions
defined in the system sys return a status other than
GSL_SUCCESS
the step will be aborted. In that case, the
elements of y will be restored to their prestep values and the
error code from the usersupplied function will be returned. Failure
may be due to a singularity in the system or too large stepsize
h. In that case the step should be attempted again with a
smaller stepsize, e.g. h/2.
If the driver object is not appropriately set via
gsl_odeiv2_step_set_driver
for those steppers that need it, the
stepping function returns GSL_EFAULT
. If the usersupplied
functions defined in the system sys returns GSL_EBADFUNC
,
the function returns immediately with the same return code. In this
case the user must call gsl_odeiv2_step_reset
before calling
this function again.
The following algorithms are available,
 Step Type: gsl_odeiv2_step_rk2

Explicit embedded RungeKutta (2, 3) method.
 Step Type: gsl_odeiv2_step_rk4

Explicit 4th order (classical) RungeKutta. Error estimation is
carried out by the step doubling method. For more efficient estimate
of the error, use the embedded methods described below.
 Step Type: gsl_odeiv2_step_rkf45

Explicit embedded RungeKuttaFehlberg (4, 5) method. This method is
a good generalpurpose integrator.
 Step Type: gsl_odeiv2_step_rkck

Explicit embedded RungeKutta CashKarp (4, 5) method.
 Step Type: gsl_odeiv2_step_rk8pd

Explicit embedded RungeKutta PrinceDormand (8, 9) method.
 Step Type: gsl_odeiv2_step_rk1imp

Implicit Gaussian first order RungeKutta. Also known as implicit
Euler or backward Euler method. Error estimation is carried out by the
step doubling method. This algorithm requires the Jacobian and
access to the driver object via
gsl_odeiv2_step_set_driver
.
 Step Type: gsl_odeiv2_step_rk2imp

Implicit Gaussian second order RungeKutta. Also known as implicit
midpoint rule. Error estimation is carried out by the step doubling
method. This stepper requires the Jacobian and access to the driver
object via
gsl_odeiv2_step_set_driver
.
 Step Type: gsl_odeiv2_step_rk4imp

Implicit Gaussian 4th order RungeKutta. Error estimation is carried
out by the step doubling method. This algorithm requires the Jacobian
and access to the driver object via
gsl_odeiv2_step_set_driver
.
 Step Type: gsl_odeiv2_step_bsimp

Implicit BulirschStoer method of Bader and Deuflhard. The method is
generally suitable for stiff problems. This stepper requires the
Jacobian.
 Step Type: gsl_odeiv2_step_msadams

A variablecoefficient linear multistep Adams method in Nordsieck
form. This stepper uses explicit AdamsBashforth (predictor) and
implicit AdamsMoulton (corrector) methods in P(EC)^m
functional iteration mode. Method order varies dynamically between 1
and 12. This stepper requires the access to the driver object via
gsl_odeiv2_step_set_driver
.
 Step Type: gsl_odeiv2_step_msbdf

A variablecoefficient linear multistep backward differentiation
formula (BDF) method in Nordsieck form. This stepper uses the explicit
BDF formula as predictor and implicit BDF formula as corrector. A
modified Newton iteration method is used to solve the system of
nonlinear equations. Method order varies dynamically between 1 and
5. The method is generally suitable for stiff problems. This stepper
requires the Jacobian and the access to the driver object via
gsl_odeiv2_step_set_driver
.
The control function examines the proposed change to the solution
produced by a stepping function and attempts to determine the optimal
stepsize for a userspecified level of error.
 Function: gsl_odeiv2_control * gsl_odeiv2_control_standard_new (double eps_abs, double eps_rel, double a_y, double a_dydt)

The standard control object is a four parameter heuristic based on
absolute and relative errors eps_abs and eps_rel, and
scaling factors a_y and a_dydt for the system state
y(t) and derivatives y'(t) respectively.
The stepsize adjustment procedure for this method begins by computing
the desired error level D_i for each component,
D_i = eps_abs + eps_rel * (a_y y_i + a_dydt h y\prime_i)
and comparing it with the observed error E_i = yerr_i. If the
observed error E exceeds the desired error level D by more
than 10% for any component then the method reduces the stepsize by an
appropriate factor,
h_new = h_old * S * (E/D)^(1/q)
where q is the consistency order of the method (e.g. q=4 for
4(5) embedded RK), and S is a safety factor of 0.9. The ratio
E/D is taken to be the maximum of the ratios
E_i/D_i.
If the observed error E is less than 50% of the desired error
level D for the maximum ratio E_i/D_i then the algorithm
takes the opportunity to increase the stepsize to bring the error in
line with the desired level,
h_new = h_old * S * (E/D)^(1/(q+1))
This encompasses all the standard error scaling methods. To avoid
uncontrolled changes in the stepsize, the overall scaling factor is
limited to the range 1/5 to 5.
 Function: gsl_odeiv2_control * gsl_odeiv2_control_y_new (double eps_abs, double eps_rel)

This function creates a new control object which will keep the local
error on each step within an absolute error of eps_abs and
relative error of eps_rel with respect to the solution y_i(t).
This is equivalent to the standard control object with a_y=1 and
a_dydt=0.
 Function: gsl_odeiv2_control * gsl_odeiv2_control_yp_new (double eps_abs, double eps_rel)

This function creates a new control object which will keep the local
error on each step within an absolute error of eps_abs and
relative error of eps_rel with respect to the derivatives of the
solution y'_i(t). This is equivalent to the standard control
object with a_y=0 and a_dydt=1.
 Function: gsl_odeiv2_control * gsl_odeiv2_control_scaled_new (double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim)

This function creates a new control object which uses the same algorithm
as
gsl_odeiv2_control_standard_new
but with an absolute error
which is scaled for each component by the array scale_abs.
The formula for D_i for this control object is,
D_i = eps_abs * s_i + eps_rel * (a_y y_i + a_dydt h y\prime_i)
where s_i is the ith component of the array scale_abs.
The same error control heuristic is used by the Matlab ODE suite.
 Function: gsl_odeiv2_control * gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type * T)

This function returns a pointer to a newly allocated instance of a
control function of type T. This function is only needed for
defining new types of control functions. For most purposes the standard
control functions described above should be sufficient.
 Function: int gsl_odeiv2_control_init (gsl_odeiv2_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt)

This function initializes the control function c with the
parameters eps_abs (absolute error), eps_rel (relative
error), a_y (scaling factor for y) and a_dydt (scaling
factor for derivatives).
 Function: void gsl_odeiv2_control_free (gsl_odeiv2_control * c)

This function frees all the memory associated with the control function
c.
 Function: int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * c, gsl_odeiv2_step * s, const double y[], const double yerr[], const double dydt[], double * h)

This function adjusts the stepsize h using the control function
c, and the current values of y, yerr and dydt.
The stepping function step is also needed to determine the order
of the method. If the error in the yvalues yerr is found to be
too large then the stepsize h is reduced and the function returns
GSL_ODEIV_HADJ_DEC
. If the error is sufficiently small then
h may be increased and GSL_ODEIV_HADJ_INC
is returned. The
function returns GSL_ODEIV_HADJ_NIL
if the stepsize is
unchanged. The goal of the function is to estimate the largest
stepsize which satisfies the userspecified accuracy requirements for
the current point.
 Function: const char * gsl_odeiv2_control_name (const gsl_odeiv2_control * c)

This function returns a pointer to the name of the control function.
For example,
printf ("control method is '%s'\n",
gsl_odeiv2_control_name (c));
would print something like control method is 'standard'
 Function: int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * c, const double y, const double dydt, const double h, const size_t ind, double * errlev)

This function calculates the desired error level of the indth component to errlev. It requires the value (y) and value of the derivative (dydt) of the component, and the current step size h.
 Function: int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * c, const gsl_odeiv2_driver * d)

This function sets a pointer of the driver object d for control
object c.
The evolution function combines the results of a stepping function and
control function to reliably advance the solution forward one step
using an acceptable stepsize.
 Function: gsl_odeiv2_evolve * gsl_odeiv2_evolve_alloc (size_t dim)

This function returns a pointer to a newly allocated instance of an
evolution function for a system of dim dimensions.
 Function: int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con, gsl_odeiv2_step * step, const gsl_odeiv2_system * sys, double * t, double t1, double * h, double y[])

This function advances the system (e, sys) from time
t and position y using the stepping function step.
The new time and position are stored in t and y on output.
The initial stepsize is taken as h. The control function
con is applied to check whether the local error estimated by the
stepping function step using stepsize h exceeds the
required error tolerance. If the error is too high, the step is
retried by calling step with a decreased stepsize. This process
is continued until an acceptable stepsize is found. An estimate of
the local error for the step can be obtained from the components of
the array e>yerr[]
.
If the usersupplied functions defined in the system sys returns
GSL_EBADFUNC
, the function returns immediately with the same
return code. In this case the user must call
gsl_odeiv2_step_reset
and
gsl_odeiv2_evolve_reset
before calling this function again.
Otherwise, if the usersupplied functions defined in the system
sys or the stepping function step return a status other
than GSL_SUCCESS
, the step is retried with a decreased
stepsize. If the stepsize decreases below machine precision, a
status of GSL_FAILURE
is returned if the user functions
returned GSL_SUCCESS
. Otherwise the value returned by user
function is returned. If no acceptable step can be made, t and
y will be restored to their prestep values and h contains
the final attempted stepsize.
If the step is successful the function returns a suggested stepsize
for the next step in h. The maximum time t1 is guaranteed
not to be exceeded by the timestep. On the final timestep the value
of t will be set to t1 exactly.
 Function: int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con, gsl_odeiv2_step * step, const gsl_odeiv2_system * sys, double * t, const double h, double y[])

This function advances the ODEsystem (e, sys, con)
from time t and position y using the stepping function
step by a specified step size h. If the local error
estimated by the stepping function exceeds the desired error level,
the step is not taken and the function returns
GSL_FAILURE
. Otherwise the value returned by user function is
returned.
 Function: int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e)

This function resets the evolution function e. It should be used
whenever the next use of e will not be a continuation of a
previous step.
 Function: void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e)

This function frees all the memory associated with the evolution function
e.
 Function: int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e, const gsl_odeiv2_driver * d)

This function sets a pointer of the driver object d for evolve
object e.
If a system has discontinuous changes in the derivatives at known
points, it is advisable to evolve the system between each discontinuity
in sequence. For example, if a stepchange in an external driving
force occurs at times t_a, t_b and t_c then evolution
should be carried out over the ranges (t_0,t_a),
(t_a,t_b), (t_b,t_c), and (t_c,t_1) separately
and not directly over the range (t_0,t_1).
The driver object is a high level wrapper that combines the evolution,
control and stepper objects for easy use.
 Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel)

 Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel)

 Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel, const double a_y, const double a_dydt)

 Function: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel, const double a_y, const double a_dydt, const double scale_abs[])

These functions return a pointer to a newly allocated instance of a
driver object. The functions automatically allocate and initialise the
evolve, control and stepper objects for ODE system sys using
stepper type T. The initial step size is given in
hstart. The rest of the arguments follow the syntax and
semantics of the control functions with same name
(
gsl_odeiv2_control_*_new
).
 Function: int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin)

The function sets a minimum for allowed step size hmin for
driver d. Default value is 0.
 Function: int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax)

The function sets a maximum for allowed step size hmax for
driver d. Default value is
GSL_DBL_MAX
.
 Function: int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d, const unsigned long int nmax)

The function sets a maximum for allowed number of steps nmax for
driver d. Default value of 0 sets no limit for steps.
 Function: int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double * t, const double t1, double y[])

This function evolves the driver system d from t to
t1. Initially vector y should contain the values of
dependent variables at point t. If the function is unable to
complete the calculation, an error code from
gsl_odeiv2_evolve_apply
is returned, and t and y
contain the values from last successful step.
If maximum number of steps is reached, a value of GSL_EMAXITER
is returned. If the step size drops below minimum value, the function
returns with GSL_ENOPROG
. If the usersupplied functions
defined in the system sys returns GSL_EBADFUNC
, the
function returns immediately with the same return code. In this case
the user must call gsl_odeiv2_driver_reset
before calling this
function again.
 Function: int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double * t, const double h, const unsigned long int n, double y[])

This function evolves the driver system d from t with
n steps of size h. If the function is unable to complete
the calculation, an error code from
gsl_odeiv2_evolve_apply_fixed_step
is returned, and t and
y contain the values from last successful step.
 Function: int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d)

This function resets the evolution and stepper objects.
 Function: int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart)

The routine resets the evolution and stepper objects and sets new
initial step size to hstart. This function can be used e.g. to
change the direction of integration.
 Function: int gsl_odeiv2_driver_free (gsl_odeiv2_driver * d)

This function frees the driver object, and the related evolution,
stepper and control objects.
The following program solves the secondorder nonlinear Van der Pol
oscillator equation,
u"(t) + \mu u'(t) (u(t)^2  1) + u(t) = 0
This can be converted into a first order system suitable for use with
the routines described in this chapter by introducing a separate
variable for the velocity, v = u'(t),
u' = v
v' = u + \mu v (1u^2)
The program begins by defining functions for these derivatives and
their Jacobian. The main function uses driver level functions to solve
the problem. The program evolves the solution from (u, v) = (1,
0) at t=0 to t=100. The stepsize h is
automatically adjusted by the controller to maintain an absolute
accuracy of
10^{6} in the function values (u, v).
The loop in the example prints the solution at the points
t_i = 1, 2, \dots, 100.
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
int
func (double t, const double y[], double f[],
void *params)
{
double mu = *(double *)params;
f[0] = y[1];
f[1] = y[0]  mu*y[1]*(y[0]*y[0]  1);
return GSL_SUCCESS;
}
int
jac (double t, const double y[], double *dfdy,
double dfdt[], void *params)
{
double mu = *(double *)params;
gsl_matrix_view dfdy_mat
= gsl_matrix_view_array (dfdy, 2, 2);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, 0.0);
gsl_matrix_set (m, 0, 1, 1.0);
gsl_matrix_set (m, 1, 0, 2.0*mu*y[0]*y[1]  1.0);
gsl_matrix_set (m, 1, 1, mu*(y[0]*y[0]  1.0));
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
int
main (void)
{
double mu = 10;
gsl_odeiv2_system sys = {func, jac, 2, &mu};
gsl_odeiv2_driver * d =
gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
1e6, 1e6, 0.0);
int i;
double t = 0.0, t1 = 100.0;
double y[2] = { 1.0, 0.0 };
for (i = 1; i <= 100; i++)
{
double ti = i * t1 / 100.0;
int status = gsl_odeiv2_driver_apply (d, &t, ti, y);
if (status != GSL_SUCCESS)
{
printf ("error, return value=%d\n", status);
break;
}
printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
}
gsl_odeiv2_driver_free (d);
return 0;
}
The user can work with the lower level functions directly, as in
the following example. In this case an intermediate result is printed
after each successful step instead of equidistant time points.
int
main (void)
{
const gsl_odeiv2_step_type * T
= gsl_odeiv2_step_rk8pd;
gsl_odeiv2_step * s
= gsl_odeiv2_step_alloc (T, 2);
gsl_odeiv2_control * c
= gsl_odeiv2_control_y_new (1e6, 0.0);
gsl_odeiv2_evolve * e
= gsl_odeiv2_evolve_alloc (2);
double mu = 10;
gsl_odeiv2_system sys = {func, jac, 2, &mu};
double t = 0.0, t1 = 100.0;
double h = 1e6;
double y[2] = { 1.0, 0.0 };
while (t < t1)
{
int status = gsl_odeiv2_evolve_apply (e, c, s,
&sys,
&t, t1,
&h, y);
if (status != GSL_SUCCESS)
break;
printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
}
gsl_odeiv2_evolve_free (e);
gsl_odeiv2_control_free (c);
gsl_odeiv2_step_free (s);
return 0;
}
For functions with multiple parameters, the appropriate information
can be passed in through the params argument in
gsl_odeiv2_system
definition (mu in this example) by using
a pointer to a struct.
It is also possible to work with a nonadaptive integrator, using only
the stepping function itself,
gsl_odeiv2_driver_apply_fixed_step
or
gsl_odeiv2_evolve_apply_fixed_step
. The following program uses
the driver level function, with fourthorder
RungeKutta stepping function with a fixed stepsize of
0.001.
int
main (void)
{
double mu = 10;
gsl_odeiv2_system sys = { func, jac, 2, &mu };
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4,
1e3, 1e8, 1e8);
double t = 0.0;
double y[2] = { 1.0, 0.0 };
int i, s;
for (i = 0; i < 100; i++)
{
s = gsl_odeiv2_driver_apply_fixed_step (d, &t, 1e3, 1000, y);
if (s != GSL_SUCCESS)
{
printf ("error: driver returned %d\n", s);
break;
}
printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
}
gsl_odeiv2_driver_free (d);
return s;
}

Ascher, U.M., Petzold, L.R., Computer Methods for Ordinary
Differential and DifferentialAlgebraic Equations, SIAM,
Philadelphia, 1998.

Hairer, E., Norsett, S. P., Wanner, G., Solving Ordinary Differential
Equations I: Nonstiff Problems, Springer, Berlin, 1993.

Hairer, E., Wanner, G., Solving Ordinary Differential
Equations II: Stiff and DifferentialAlgebraic Problems,
Springer, Berlin, 1996.
Many of the basic RungeKutta formulas can be found in the Handbook of
Mathematical Functions,

Abramowitz & Stegun (eds.), Handbook of Mathematical Functions,
Section 25.5.
The implicit BulirschStoer algorithm bsimp
is described in the
following paper,

G. Bader and P. Deuflhard, "A SemiImplicit MidPoint Rule for Stiff
Systems of Ordinary Differential Equations.", Numer. Math. 41, 373398,
1983.
The Adams and BDF multistep methods msadams
and msbdf
are based on the following articles,

G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the
Numerical Solution of Ordinary Differential Equations.",
ACM Trans. Math. Software, 1, 7196, 1975.

P. N. Brown, G. D. Byrne and A. C. Hindmarsh, "VODE: A
Variablecoefficient ODE Solver.", SIAM J. Sci. Stat. Comput. 10,
10381051, 1989.

A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban,
D. E. Shumaker and C. S. Woodward, "SUNDIALS: Suite of
Nonlinear and Differential/Algebraic Equation Solvers.", ACM
Trans. Math. Software 31, 363396, 2005.
Go to the first, previous, next, last section, table of contents.