ode-initval.texi

开放gsl矩阵运算

@cindex differential equations, initial value problems
@cindex initial value problems, differential equations
@cindex ordinary differential equations, initial value problem
@cindex ODEs, initial value problems
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 @file{gsl_odeiv.h}.

@menu
* Defining the ODE System::     
* Stepping Functions::          
* Adaptive Step-size Control::  
* Evolution::                   
* ODE Example programs::        
* ODE References and Further Reading::  
@end menu

@node Defining the ODE System
@section Defining the ODE System

The routines solve the general @math{n}-dimensional first-order system,

@tex
\beforedisplay
$$
{dy_i(t) \over dt} = f_i (t, y_1(t), \dots y_n(t))
$$
\afterdisplay
@end tex
@ifinfo
@example
dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t))
@end example
@end ifinfo
@noindent
for @math{i = 1, \dots, n}.  The stepping functions rely on the vector
of derivatives @math{f_i} and the Jacobian matrix, 
@c{$J_{ij} = \partial f_i(t, y(t)) / \partial y_j$}
@math{J_@{ij@} = df_i(t,y(t)) / dy_j}. 
A system of equations is defined using the @code{gsl_odeiv_system}
datatype.

@deftp {Data Type} gsl_odeiv_system
This data type defines a general ODE system with arbitrary parameters.

@table @code
@item int (* @var{function}) (double t, const double y[], double dydt[], void * params)
This function should store the elements of
@c{$f(t,y,\hbox{\it params})$}
@math{f(t,y,params)} in the array @var{dydt},
for arguments (@var{t},@var{y}) and parameters @var{params}

@item  int (* @var{jacobian}) (double t, const double y[], double * dfdy, double dfdt[], void * params);
This function should store the elements of @math{f(t,y,params)} in
the array @var{dfdt} and the Jacobian matrix @c{$J_{ij}$}
@math{J_@{ij@}} in the the array
@var{dfdy} regarded as a row-ordered matrix @code{J(i,j) = dfdy[i * dim + j]}
where @code{dim} is the dimension of the system.

@item size_t @var{dimension};
This is the dimension of the system of equations

@item void * @var{params}
This is a pointer to the arbitrary parameters of the system.
@end table
@end deftp

@node Stepping Functions
@section Stepping Functions

The lowest level components are the @dfn{stepping functions} which
advance a solution from time @math{t} to @math{t+h} for a fixed
step-size @math{h} and estimate the resulting local error.

@deftypefun {gsl_odeiv_step *} gsl_odeiv_step_alloc (const gsl_odeiv_step_type * @var{T}, size_t @var{dim})
This function returns a pointer to a newly allocated instance of a
stepping function of type @var{T} for a system of @var{dim} dimensions.
@end deftypefun

@deftypefun int gsl_odeiv_step_reset (gsl_odeiv_step * @var{s})
This function resets the stepping function @var{s}.  It should be used
whenever the next use of @var{s} will not be a continuation of a
previous step.
@end deftypefun

@deftypefun void gsl_odeiv_step_free (gsl_odeiv_step * @var{s})
This function frees all the memory associated with the stepping function
@var{s}.
@end deftypefun

@deftypefun {const char *} gsl_odeiv_step_name (const gsl_odeiv_step * @var{s})
This function returns a pointer to the name of the stepping function.
For example,

@example
printf("step method is '%s'\n",
        gsl_odeiv_step_name (s));
@end example

@noindent
would print something like @code{step method is 'rk4'}.
@end deftypefun

@deftypefun {unsigned int} gsl_odeiv_step_order (const gsl_odeiv_step * @var{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.
@end deftypefun

@deftypefun int gsl_odeiv_step_apply (gsl_odeiv_step * @var{s}, double @var{t}, double @var{h}, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * @var{dydt})
This function applies the stepping function @var{s} to the system of
equations defined by @var{dydt}, using the step size @var{h} to advance
the system from time @var{t} and state @var{y} to time @var{t}+@var{h}.
The new state of the system is stored in @var{y} on output, with an
estimate of the absolute error in each component stored in @var{yerr}.
If the argument @var{dydt_in} is not null it should point an array
containing the derivatives for the system at time @var{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 @var{t}+@var{h} will
be stored in @var{dydt_out} if it is not null.
@end deftypefun

The following algorithms are available,

@deffn {Step Type} gsl_odeiv_step_rk2
@cindex RK2, Runge-Kutta Method
@cindex Runge-Kutta methods, ordinary differential equations
Embedded 2nd order Runge-Kutta with 3rd order error estimate.
@end deffn

@deffn {Step Type} gsl_odeiv_step_rk4
@cindex RK4, Runge-Kutta Method
4th order (classical) Runge-Kutta.
@end deffn

@deffn {Step Type} gsl_odeiv_step_rkf45
@cindex Fehlberg method, differential equations
@cindex RKF45, Runge-Kutta-Fehlberg method
Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error
estimate.  This method is a good general-purpose integrator.
@end deffn

@deffn {Step Type} gsl_odeiv_step_rkck 
@cindex Runge-Kutta Cash-Karp method
@cindex Cash-Karp, Runge-Kutta method
Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error
estimate.
@end deffn

@deffn {Step Type} gsl_odeiv_step_rk8pd  
@cindex Runge-Kutta Prince-Dormand method
@cindex Prince-Dormand, Runge-Kutta method
Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order
error estimate.
@end deffn

@deffn {Step Type} gsl_odeiv_step_rk2imp  
Implicit 2nd order Runge-Kutta at Gaussian points
@end deffn

@deffn {Step Type} gsl_odeiv_step_rk4imp  
Implicit 4th order Runge-Kutta at Gaussian points
@end deffn

@deffn {Step Type} gsl_odeiv_step_bsimp   
@cindex Bulirsch-Stoer method
@cindex Bader and Deuflhard, Bulirsch-Stoer method.
@cindex Deuflhard and Bader, Bulirsch-Stoer method.
Implicit Bulirsch-Stoer method of Bader and Deuflhard.
@end deffn

@deffn {Step Type} gsl_odeiv_step_gear1 
@cindex Gear method, differential equations
M=1 implicit Gear method 
@end deffn

@deffn {Step Type} gsl_odeiv_step_gear2 
M=2 implicit Gear method
@end deffn

@node Adaptive Step-size Control
@section Adaptive Step-size Control
@cindex Adaptive step-size control, differential equations

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.

@deftypefun {gsl_odeiv_control *} gsl_odeiv_control_standard_new (double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt})
The standard control object is a four parameter heuristic based on
absolute and relative errors @var{eps_abs} and @var{eps_rel}, and
scaling factors @var{a_y} and @var{a_dydt} for the system state
@math{y(t)} and derivatives @math{y'(t)} respectively.

The step-size adjustment procedure for this method begins by computing
the desired error level @math{D_i} for each component,

@tex
\beforedisplay
$$
D_i = \epsilon_{abs} + \epsilon_{rel} * (a_{y} |y_i| + a_{dydt} h |y'_i|)
$$
\afterdisplay
@end tex
@ifinfo
@example
D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)
@end example
@end ifinfo
@noindent
and comparing it with the observed error @math{E_i = |yerr_i|}.  If the
observed error @var{E} exceeds the desired error level @var{D} by more
than 10% for any component then the method reduces the step-size by an
appropriate factor,

@tex
\beforedisplay
$$
h_{new} = h_{old} * S * (D/E)^{1/q}
$$
\afterdisplay
@end tex
@ifinfo
@example
h_new = h_old * S * (D/E)^(1/q)
@end example
@end ifinfo
@noindent
where @math{q} is the consistency order of method (e.g. @math{q=4} for
4(5) embedded RK), and @math{S} is a safety factor of 0.9. The ratio
@math{D/E} is taken to be the maximum of the ratios
@math{D_i/E_i}. 

If the observed error @math{E} is less than 50% of the desired error
level @var{D} for the maximum ratio @math{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,

@tex
\beforedisplay
$$
h_{new} = h_{old} * S * (E/D)^{1/(q+1)}
$$
\afterdisplay
@end tex
@ifinfo
@example
h_new = h_old * S * (E/D)^(1/(q+1))
@end example
@end ifinfo
@noindent
This encompasses all the standard error scaling methods.
@end deftypefun

@deftypefun {gsl_odeiv_control *} gsl_odeiv_control_y_new (double @var{eps_abs}, double @var{eps_rel})
This function creates a new control object which will keep the local
error on each step within an absolute error of @var{eps_abs} and
relative error of @var{eps_rel} with respect to the solution @math{y_i(t)}.
This is equivalent to the standard control object with @var{a_y}=1 and
@var{a_dydt}=0.
@end deftypefun

@deftypefun {gsl_odeiv_control *} gsl_odeiv_control_yp_new (double @var{eps_abs}, double @var{eps_rel})
This function creates a new control object which will keep the local
error on each step within an absolute error of @var{eps_abs} and