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 low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines, and higher-level components for adaptive step-size control. The components can be combined by the user to achieve the desired solution, with full access to any intermediate steps.

These functions are declared in the header file ``gsl_odeiv.h'`.

The routines solve the general n-dimensional first-order 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_odeiv_system`

datatype.

__Data Type:__**gsl_odeiv_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` `int (*`

`jacobian`) (double t, const double y[], double * dfdy, double dfdt[], void * params);-
This function should store the vector elements
df_i(t,y,params)/dt in the array
`dfdt`and the Jacobian matrix J_{ij} in the the array`dfdy`regarded as a row-ordered matrix`J(i,j) = dfdy[i * dim + j]`

where`dim`

is the dimension of the system. Some of the simpler solver algorithms do not make use of the Jacobian matrix, so it is not always strictly necessary to provide it (this element of the struct can be replace by a null pointer). However, it is useful to provide the Jacobian to allow the solver algorithms to be interchanged -- the best algorithms make use of the Jacobian. `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
step-size h and estimate the resulting local error.

__Function:__gsl_odeiv_step ***gsl_odeiv_step_alloc***(const gsl_odeiv_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.

__Function:__int**gsl_odeiv_step_reset***(gsl_odeiv_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_odeiv_step_free***(gsl_odeiv_step **`s`)-
This function frees all the memory associated with the stepping function
`s`.

__Function:__const char ***gsl_odeiv_step_name***(const gsl_odeiv_step **`s`)-
This function returns a pointer to the name of the stepping function.
For example,
printf("step method is '%s'\n", gsl_odeiv_step_name (s));

would print something like

`step method is 'rk4'`

.

__Function:__unsigned int**gsl_odeiv_step_order***(const gsl_odeiv_step **`s`)- This function returns the order of the stepping function on the previous step. This order can vary if the stepping function itself is adaptive.

__Function:__int**gsl_odeiv_step_apply***(gsl_odeiv_step **`s`, double`t`, double`h`, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system *`dydt`)-
This function applies the stepping function
`s`to the system of equations defined by`dydt`, using the step size`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 following algorithms are available,

__Step Type:__**gsl_odeiv_step_rk2**- Embedded 2nd order Runge-Kutta with 3rd order error estimate.

__Step Type:__**gsl_odeiv_step_rk4**- 4th order (classical) Runge-Kutta.

__Step Type:__**gsl_odeiv_step_rkf45**- Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error estimate. This method is a good general-purpose integrator.

__Step Type:__**gsl_odeiv_step_rkck**- Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error estimate.

__Step Type:__**gsl_odeiv_step_rk8pd**- Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimate.

__Step Type:__**gsl_odeiv_step_rk2imp**- Implicit 2nd order Runge-Kutta at Gaussian points

__Step Type:__**gsl_odeiv_step_rk4imp**- Implicit 4th order Runge-Kutta at Gaussian points

__Step Type:__**gsl_odeiv_step_bsimp**- Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.

__Step Type:__**gsl_odeiv_step_gear1**- M=1 implicit Gear method

__Step Type:__**gsl_odeiv_step_gear2**- M=2 implicit Gear method

The control function examines the proposed change to the solution and its error estimate produced by a stepping function and attempts to determine the optimal step-size for a user-specified level of error.

__Function:__gsl_odeiv_control ***gsl_odeiv_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 step-size 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'_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 step-size by an appropriate factor,h_new = h_old * S * (D/E)^(1/q)

where q is the consistency order of method (e.g. q=4 for 4(5) embedded RK), and S is a safety factor of 0.9. The ratio D/E is taken to be the maximum of the ratios D_i/E_i.

If the observed error E is less than 50% of the desired error level

`D`for the maximum ratio D_i/E_i then the algorithm takes the opportunity to increase the step-size 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.

__Function:__gsl_odeiv_control ***gsl_odeiv_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_odeiv_control ***gsl_odeiv_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_odeiv_control ***gsl_odeiv_control_alloc***(const gsl_odeiv_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_odeiv_control_init***(gsl_odeiv_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_odeiv_control_free***(gsl_odeiv_control **`c`)-
This function frees all the memory associated with the control function
`c`.

__Function:__int**gsl_odeiv_control_hadjust***(gsl_odeiv_control **`c`, gsl_odeiv_step *`s`, const double y0[], const double yerr[], const double dydt[], double *`h`)-
This function adjusts the step-size
`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 y-values`yerr`is found to be too large then the step-size`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 step-size is unchanged. The goal of the function is to estimate the largest step-size which satisfies the user-specified accuracy requirements for the current point.

__Function:__const char ***gsl_odeiv_control_name***(const gsl_odeiv_control **`c`)-
This function returns a pointer to the name of the control function.
For example,
printf("control method is '%s'\n", gsl_odeiv_control_name (c));

would print something like

`control method is 'standard'`

The highest level of the system is the evolution function which combines the results of a stepping function and control function to reliably advance the solution forward over an interval (t_0, t_1). If the control function signals that the step-size should be decreased the evolution function backs out of the current step and tries the proposed smaller step-size. This is process is continued until an acceptable step-size is found.

__Function:__gsl_odeiv_evolve ***gsl_odeiv_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_odeiv_evolve_apply***(gsl_odeiv_evolve **`e`, gsl_odeiv_control *`con`, gsl_odeiv_step *`step`, const gsl_odeiv_system *`dydt`, double *`t`, double`t1`, double *`h`, double y[])-
This function advances the system (
`e`,`dydt`) 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 step-size is taken as`h`, but this will be modified using the control function`c`to achieve the appropriate error bound if necessary. The routine may make several calls to`step`in order to determine the optimum step-size. If the step-size has been changed the value of`h`will be modified on output. The maximum time`t1`is guaranteed not to be exceeded by the time-step. On the final time-step the value of`t`will be set to`t1`exactly.

__Function:__int**gsl_odeiv_evolve_reset***(gsl_odeiv_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_odeiv_evolve_free***(gsl_odeiv_evolve **`e`)-
This function frees all the memory associated with the evolution function
`e`.

The following program solves the second-order nonlinear Van der Pol oscillator equation,

x"(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0

This can be converted into a first order system suitable for use with the library by introducing a separate variable for the velocity, y = x'(t),

x' = y y' = -x + \mu y (1-x^2)

The program begins by defining functions for these derivatives and their Jacobian,

#include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_odeiv.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) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0); gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2); double mu = 10; gsl_odeiv_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-6; double y[2] = { 1.0, 0.0 }; while (t < t1) { int status = gsl_odeiv_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_odeiv_evolve_free(e); gsl_odeiv_control_free(c); gsl_odeiv_step_free(s); return 0; }

The main loop of the program evolves the solution from (y, y') =
(1, 0) at t=0 to t=100. The step-size h is
automatically adjusted by the controller to maintain an absolute
accuracy of
10^{-6} in the function values `y`.

To obtain the values at regular intervals, rather than the variable spacings chosen by the control function, the main loop can be modified to advance the solution from one point to the next. For example, the following main loop prints the solution at the fixed points t = 0, 1, 2, \dots, 100,

for (i = 1; i <= 100; i++) { double ti = i * t1 / 100.0; while (t < ti) { gsl_odeiv_evolve_apply (e, c, s, &sys, &t, ti, &h, y); } printf("%.5e %.5e %.5e\n", t, y[0], y[1]); }

It is also possible to work with a non-adaptive integrator, using only
the stepping function itself. The following program uses the `rk4`

fourth-order Runge-Kutta stepping function with a fixed stepsize of
0.01,

int main (void) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); double mu = 10; gsl_odeiv_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-2; double y[2] = { 1.0, 0.0 }, y_err[2]; double dfdy[4], dydt_in[2], dydt_out[2]; /* initialise dydt_in */ GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in); while (t < t1) { int status = gsl_odeiv_step_apply (s, t, h, y, y_err, dydt_in, dydt_out, &sys); if (status != GSL_SUCCESS) break; dydt_in[0] = dydt_out[0]; dydt_in[1] = dydt_out[1]; t += h; printf("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv_step_free(s); return 0; }

The derivatives and jacobian must be initialised for the starting point
t=0 before the first step is taken. Subsequent steps use the
output derivatives `dydt_out` as inputs to the next step by copying
their values into `dydt_in`.

Many of the the basic Runge-Kutta formulas can be found in the Handbook of Mathematical Functions,

- Abramowitz & Stegun (eds.), Handbook of Mathematical Functions, Section 25.5.

The implicit Bulirsch-Stoer algorithm `bsimp`

is described in the
following paper,

- G. Bader and P. Deuflhard, "A Semi-Implicit Mid-Point Rule for Stiff Systems of Ordinary Differential Equations.", Numer. Math. 41, 373-398, 1983.

Go to the first, previous, next, last section, table of contents.