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.
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.
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.
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.
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,
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.
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.
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.
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.
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.
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.
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.
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.