src/parsec/disk-image/parsec/parsec-benchmark/pkgs/libs/gsl/src/doc/ode-initval.texi - public/gem5-resources - Git at Google

 @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 (* function) (double t, const double y[], double dydt[], void * params)
 This function should store the vector elements
 @c{$f_i(t,y,\hbox{\it params})$}
 @math{f_i(t,y,params)} in the array @var{dydt},
 for arguments (@var{t},@var{y}) and parameters @var{params}.
 The function should return @code{GSL_SUCCESS} if the calculation
 was completed successfully. Any other return value indicates
 an error.

 @item  int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
 This function should store the vector of derivative elements
 @c{$\partial f_i(t,y,params) / \partial t$}
 @math{df_i(t,y,params)/dt} in the array @var{dfdt} and the
 Jacobian matrix @c{$J_{ij}$}
 @math{J_@{ij@}} in the array
 @var{dfdy}, regarded as a row-ordered matrix @code{J(i,j) = dfdy[i * dimension + j]}
 where @code{dimension} is the dimension of the system.
 The function should return @code{GSL_SUCCESS} if the calculation
 was completed successfully.  Any other return value indicates
 an error.

 Some of the simpler solver algorithms do not make use of the Jacobian
 matrix, so it is not always strictly necessary to provide it (the
 @code{jacobian} element of the struct can be replaced by a null pointer
 for those algorithms).  However, it is useful to provide the Jacobian to allow
 the solver algorithms to be interchanged---the best algorithms make use
 of the Jacobian.

 @item size_t dimension;
 This is the dimension of the system of equations.

 @item void * 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 @var{y}[], double @var{yerr}[], const double @var{dydt_in}[], double @var{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.

 If the user-supplied functions defined in the system @var{dydt} return a
 status other than @code{GSL_SUCCESS} the step will be aborted.  In this
 case, the elements of @var{y} will be restored to their pre-step values
 and the error code from the user-supplied function will be returned.
 The step-size @var{h} will be set to the step-size which caused the error.
 If the function is called again with a smaller step-size, e.g. @math{@var{h}/10},
 it should be possible to get closer to any singularity.
 To distinguish between error codes from the user-supplied functions and
 those from @code{gsl_odeiv_step_apply} itself, any user-defined return
 values should be distinct from the standard GSL error codes.
 @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 Runge-Kutta (2, 3) method.
 @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 Runge-Kutta-Fehlberg (4, 5) method.  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 Runge-Kutta Cash-Karp (4, 5) method.
 @end deffn

 @deffn {Step Type} gsl_odeiv_step_rk8pd
 @cindex Runge-Kutta Prince-Dormand method
 @cindex Prince-Dormand, Runge-Kutta method
 Embedded Runge-Kutta Prince-Dormand (8,9) method.
 @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.  This algorithm
 requires the Jacobian.
 @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
 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 * (E/D)^{-1/q}
 $$
 \afterdisplay
 @end tex
 @ifinfo

 @example
 h_new = h_old * S * (E/D)^(-1/q)
 @end example

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

 If the observed error @math{E} is less than 50% of the desired error
 level @var{D} for the maximum ratio @math{E_i/D_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. To avoid
 uncontrolled changes in the stepsize, the overall scaling factor is
 limited to the range @math{1/5} to 5.
 @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
 relative error of @var{eps_rel} with respect to the derivatives of the
 solution @math{y'_i(t)}.  This is equivalent to the standard control
 object with @var{a_y}=0 and @var{a_dydt}=1.
 @end deftypefun


 @deftypefun {gsl_odeiv_control *} gsl_odeiv_control_scaled_new (double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt}, const double @var{scale_abs}[], size_t @var{dim})
 This function creates a new control object which uses the same algorithm
 as @code{gsl_odeiv_control_standard_new} but with an absolute error
 which is scaled for each component by the array @var{scale_abs}.
 The formula for @math{D_i} for this control object is,
 @tex
 \beforedisplay
 $$
 D_i = \epsilon_{abs} s_i + \epsilon_{rel} * (a_{y} |y_i| + a_{dydt} h |y'_i|)
 $$
 \afterdisplay
 @end tex
 @ifinfo

 @example
 D_i = eps_abs * s_i + eps_rel * (a_y |y_i| + a_dydt h |y'_i|)
 @end example

 @end ifinfo
 @noindent
 where @math{s_i} is the @math{i}-th component of the array @var{scale_abs}.
 The same error control heuristic is used by the Matlab @sc{ode} suite.
 @end deftypefun

 @deftypefun {gsl_odeiv_control *} gsl_odeiv_control_alloc (const gsl_odeiv_control_type * @var{T})
 This function returns a pointer to a newly allocated instance of a
 control function of type @var{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.
 @end deftypefun

 @deftypefun int gsl_odeiv_control_init (gsl_odeiv_control * @var{c}, double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt})
 This function initializes the control function @var{c} with the
 parameters @var{eps_abs} (absolute error), @var{eps_rel} (relative
 error), @var{a_y} (scaling factor for y) and @var{a_dydt} (scaling
 factor for derivatives).
 @end deftypefun

 @deftypefun void gsl_odeiv_control_free (gsl_odeiv_control * @var{c})
 This function frees all the memory associated with the control function
 @var{c}.
 @end deftypefun

 @deftypefun int gsl_odeiv_control_hadjust (gsl_odeiv_control * @var{c}, gsl_odeiv_step * @var{s}, const double @var{y}[], const double @var{yerr}[], const double @var{dydt}[], double * @var{h})
 This function adjusts the step-size @var{h} using the control function
 @var{c}, and the current values of @var{y}, @var{yerr} and @var{dydt}.
 The stepping function @var{step} is also needed to determine the order
 of the method.  If the error in the y-values @var{yerr} is found to be
 too large then the step-size @var{h} is reduced and the function returns
 @code{GSL_ODEIV_HADJ_DEC}.  If the error is sufficiently small then
 @var{h} may be increased and @code{GSL_ODEIV_HADJ_INC} is returned.  The
 function returns @code{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.
 @end deftypefun

 @deftypefun {const char *} gsl_odeiv_control_name (const gsl_odeiv_control * @var{c})
 This function returns a pointer to the name of the control function.
 For example,

 @example
 printf ("control method is '%s'\n",
         gsl_odeiv_control_name (c));
 @end example

 @noindent
 would print something like @code{control method is 'standard'}
 @end deftypefun


 @node Evolution
 @section Evolution

 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 @math{(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 process is continued until an acceptable
 step-size is found.

 @deftypefun {gsl_odeiv_evolve *} gsl_odeiv_evolve_alloc (size_t @var{dim})
 This function returns a pointer to a newly allocated instance of an
 evolution function for a system of @var{dim} dimensions.
 @end deftypefun

 @deftypefun int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * @var{e}, gsl_odeiv_control * @var{con}, gsl_odeiv_step * @var{step}, const gsl_odeiv_system * @var{dydt}, double * @var{t}, double @var{t1}, double * @var{h}, double @var{y}[])
 This function advances the system (@var{e}, @var{dydt}) from time
 @var{t} and position @var{y} using the stepping function @var{step}.
 The new time and position are stored in @var{t} and @var{y} on output.
 The initial step-size is taken as @var{h}, but this will be modified
 using the control function @var{c} to achieve the appropriate error
 bound if necessary.  The routine may make several calls to @var{step} in
 order to determine the optimum step-size. An estimate of
 the local error for the step can be obtained from the components of the array @code{@var{e}->yerr[]}.
 If the step-size has been changed the value of @var{h} will be modified on output.
 The maximum time @var{t1} is guaranteed not to be exceeded by the time-step.  On the
 final time-step the value of @var{t} will be set to @var{t1} exactly.

 If the user-supplied functions defined in the system @var{dydt} return a
 status other than @code{GSL_SUCCESS} the step will be aborted.  In this
 case,  @var{t} and @var{y} will be restored to their pre-step values
 and the error code from the user-supplied function will be returned.  To
 distinguish between error codes from the user-supplied functions and
 those from @code{gsl_odeiv_evolve_apply} itself, any user-defined return
 values should be distinct from the standard GSL error codes.
 @end deftypefun

 @deftypefun int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * @var{e})
 This function resets the evolution function @var{e}.  It should be used
 whenever the next use of @var{e} will not be a continuation of a
 previous step.
 @end deftypefun

 @deftypefun void gsl_odeiv_evolve_free (gsl_odeiv_evolve * @var{e})
 This function frees all the memory associated with the evolution function
 @var{e}.
 @end deftypefun

 @node ODE Example programs
 @section Examples
 @cindex Van der Pol oscillator, example
 The following program solves the second-order nonlinear Van der Pol
 oscillator equation,
 @tex
 \beforedisplay
 $$
 x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
 $$
 \afterdisplay
 @end tex
 @ifinfo

 @example
 x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
 @end example

 @end ifinfo
 @noindent
 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, @math{y = x'(t)},
 @tex
 \beforedisplay
 $$
 \eqalign{
 x' &= y\cr
 y' &= -x + \mu y (1-x^2)
 }
 $$
 \afterdisplay
 @end tex
 @ifinfo

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

 @end ifinfo
 @noindent
 The program begins by defining functions for these derivatives and
 their Jacobian,

 @example
 @verbatiminclude examples/ode-initval.c
 @end example

 @noindent
 For functions with multiple parameters, the appropriate information
 can be passed in through the @var{params} argument using a pointer to
 a struct.

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

 @iftex
 @sp 1
 @center @image{vdp,3.4in}
 @center Numerical solution of the Van der Pol oscillator equation
 @center using Prince-Dormand 8th order Runge-Kutta.
 @end iftex

 @noindent
 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 @math{t = 0,
 1, 2, \dots, 100},

 @example
   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]);
     @}
 @end example

 @noindent
 It is also possible to work with a non-adaptive integrator, using only
 the stepping function itself.  The following program uses the @code{rk4}
 fourth-order Runge-Kutta stepping function with a fixed stepsize of
 0.01,

 @example
 @verbatiminclude examples/odefixed.c
 @end example

 @noindent
 The derivatives must be initialized for the starting point @math{t=0}
 before the first step is taken.  Subsequent steps use the output
 derivatives @var{dydt_out} as inputs to the next step by copying their
 values into @var{dydt_in}.

 @node ODE References and Further Reading
 @section References and Further Reading

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

 @itemize @asis
 @item
 Abramowitz & Stegun (eds.), @cite{Handbook of Mathematical Functions},
 Section 25.5.
 @end itemize

 @noindent
 The implicit Bulirsch-Stoer algorithm @code{bsimp} is described in the
 following paper,

 @itemize @asis
 @item
 G. Bader and P. Deuflhard, ``A Semi-Implicit Mid-Point Rule for Stiff
 Systems of Ordinary Differential Equations.'', Numer.@: Math.@: 41, 373--398,
 1983.
 @end itemize
	@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 (* function) (double t, const double y[], double dydt[], void * params)
	This function should store the vector elements
	@c{$f_i(t,y,\hbox{\it params})$}
	@math{f_i(t,y,params)} in the array @var{dydt},
	for arguments (@var{t},@var{y}) and parameters @var{params}.
	The function should return @code{GSL_SUCCESS} if the calculation
	was completed successfully. Any other return value indicates
	an error.

	@item int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
	This function should store the vector of derivative elements
	@c{$\partial f_i(t,y,params) / \partial t$}
	@math{df_i(t,y,params)/dt} in the array @var{dfdt} and the
	Jacobian matrix @c{$J_{ij}$}
	@math{J_@{ij@}} in the array
	@var{dfdy}, regarded as a row-ordered matrix @code{J(i,j) = dfdy[i * dimension + j]}
	where @code{dimension} is the dimension of the system.
	The function should return @code{GSL_SUCCESS} if the calculation
	was completed successfully. Any other return value indicates
	an error.

	Some of the simpler solver algorithms do not make use of the Jacobian
	matrix, so it is not always strictly necessary to provide it (the
	@code{jacobian} element of the struct can be replaced by a null pointer
	for those algorithms). However, it is useful to provide the Jacobian to allow
	the solver algorithms to be interchanged---the best algorithms make use
	of the Jacobian.

	@item size_t dimension;
	This is the dimension of the system of equations.

	@item void * 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 @var{y}[], double @var{yerr}[], const double @var{dydt_in}[], double @var{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.

	If the user-supplied functions defined in the system @var{dydt} return a
	status other than @code{GSL_SUCCESS} the step will be aborted. In this
	case, the elements of @var{y} will be restored to their pre-step values
	and the error code from the user-supplied function will be returned.
	The step-size @var{h} will be set to the step-size which caused the error.
	If the function is called again with a smaller step-size, e.g. @math{@var{h}/10},
	it should be possible to get closer to any singularity.
	To distinguish between error codes from the user-supplied functions and
	those from @code{gsl_odeiv_step_apply} itself, any user-defined return
	values should be distinct from the standard GSL error codes.
	@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 Runge-Kutta (2, 3) method.
	@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 Runge-Kutta-Fehlberg (4, 5) method. 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 Runge-Kutta Cash-Karp (4, 5) method.
	@end deffn

	@deffn {Step Type} gsl_odeiv_step_rk8pd
	@cindex Runge-Kutta Prince-Dormand method
	@cindex Prince-Dormand, Runge-Kutta method
	Embedded Runge-Kutta Prince-Dormand (8,9) method.
	@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. This algorithm
	requires the Jacobian.
	@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
	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 * (E/D)^{-1/q}
	$$
	\afterdisplay
	@end tex
	@ifinfo

	@example
	h_new = h_old * S * (E/D)^(-1/q)
	@end example

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

	If the observed error @math{E} is less than 50% of the desired error
	level @var{D} for the maximum ratio @math{E_i/D_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. To avoid
	uncontrolled changes in the stepsize, the overall scaling factor is
	limited to the range @math{1/5} to 5.
	@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
	relative error of @var{eps_rel} with respect to the derivatives of the
	solution @math{y'_i(t)}. This is equivalent to the standard control
	object with @var{a_y}=0 and @var{a_dydt}=1.
	@end deftypefun


	@deftypefun {gsl_odeiv_control *} gsl_odeiv_control_scaled_new (double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt}, const double @var{scale_abs}[], size_t @var{dim})
	This function creates a new control object which uses the same algorithm
	as @code{gsl_odeiv_control_standard_new} but with an absolute error
	which is scaled for each component by the array @var{scale_abs}.
	The formula for @math{D_i} for this control object is,
	@tex
	\beforedisplay
	$$
	D_i = \epsilon_{abs} s_i + \epsilon_{rel} * (a_{y} \|y_i\| + a_{dydt} h \|y'_i\|)
	$$
	\afterdisplay
	@end tex
	@ifinfo

	@example
	D_i = eps_abs * s_i + eps_rel * (a_y \|y_i\| + a_dydt h \|y'_i\|)
	@end example

	@end ifinfo
	@noindent
	where @math{s_i} is the @math{i}-th component of the array @var{scale_abs}.
	The same error control heuristic is used by the Matlab @sc{ode} suite.
	@end deftypefun

	@deftypefun {gsl_odeiv_control } gsl_odeiv_control_alloc (const gsl_odeiv_control_type @var{T})
	This function returns a pointer to a newly allocated instance of a
	control function of type @var{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.
	@end deftypefun

	@deftypefun int gsl_odeiv_control_init (gsl_odeiv_control * @var{c}, double @var{eps_abs}, double @var{eps_rel}, double @var{a_y}, double @var{a_dydt})
	This function initializes the control function @var{c} with the
	parameters @var{eps_abs} (absolute error), @var{eps_rel} (relative
	error), @var{a_y} (scaling factor for y) and @var{a_dydt} (scaling
	factor for derivatives).
	@end deftypefun

	@deftypefun void gsl_odeiv_control_free (gsl_odeiv_control * @var{c})
	This function frees all the memory associated with the control function
	@var{c}.
	@end deftypefun

	@deftypefun int gsl_odeiv_control_hadjust (gsl_odeiv_control * @var{c}, gsl_odeiv_step * @var{s}, const double @var{y}[], const double @var{yerr}[], const double @var{dydt}[], double * @var{h})
	This function adjusts the step-size @var{h} using the control function
	@var{c}, and the current values of @var{y}, @var{yerr} and @var{dydt}.
	The stepping function @var{step} is also needed to determine the order
	of the method. If the error in the y-values @var{yerr} is found to be
	too large then the step-size @var{h} is reduced and the function returns
	@code{GSL_ODEIV_HADJ_DEC}. If the error is sufficiently small then
	@var{h} may be increased and @code{GSL_ODEIV_HADJ_INC} is returned. The
	function returns @code{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.
	@end deftypefun

	@deftypefun {const char } gsl_odeiv_control_name (const gsl_odeiv_control @var{c})
	This function returns a pointer to the name of the control function.
	For example,

	@example
	printf ("control method is '%s'\n",
	gsl_odeiv_control_name (c));
	@end example

	@noindent
	would print something like @code{control method is 'standard'}
	@end deftypefun


	@node Evolution
	@section Evolution

	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 @math{(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 process is continued until an acceptable
	step-size is found.

	@deftypefun {gsl_odeiv_evolve *} gsl_odeiv_evolve_alloc (size_t @var{dim})
	This function returns a pointer to a newly allocated instance of an
	evolution function for a system of @var{dim} dimensions.
	@end deftypefun

	@deftypefun int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * @var{e}, gsl_odeiv_control * @var{con}, gsl_odeiv_step * @var{step}, const gsl_odeiv_system * @var{dydt}, double * @var{t}, double @var{t1}, double * @var{h}, double @var{y}[])
	This function advances the system (@var{e}, @var{dydt}) from time
	@var{t} and position @var{y} using the stepping function @var{step}.
	The new time and position are stored in @var{t} and @var{y} on output.
	The initial step-size is taken as @var{h}, but this will be modified
	using the control function @var{c} to achieve the appropriate error
	bound if necessary. The routine may make several calls to @var{step} in
	order to determine the optimum step-size. An estimate of
	the local error for the step can be obtained from the components of the array @code{@var{e}->yerr[]}.
	If the step-size has been changed the value of @var{h} will be modified on output.
	The maximum time @var{t1} is guaranteed not to be exceeded by the time-step. On the
	final time-step the value of @var{t} will be set to @var{t1} exactly.

	If the user-supplied functions defined in the system @var{dydt} return a
	status other than @code{GSL_SUCCESS} the step will be aborted. In this
	case, @var{t} and @var{y} will be restored to their pre-step values
	and the error code from the user-supplied function will be returned. To
	distinguish between error codes from the user-supplied functions and
	those from @code{gsl_odeiv_evolve_apply} itself, any user-defined return
	values should be distinct from the standard GSL error codes.
	@end deftypefun

	@deftypefun int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * @var{e})
	This function resets the evolution function @var{e}. It should be used
	whenever the next use of @var{e} will not be a continuation of a
	previous step.
	@end deftypefun

	@deftypefun void gsl_odeiv_evolve_free (gsl_odeiv_evolve * @var{e})
	This function frees all the memory associated with the evolution function
	@var{e}.
	@end deftypefun

	@node ODE Example programs
	@section Examples
	@cindex Van der Pol oscillator, example
	The following program solves the second-order nonlinear Van der Pol
	oscillator equation,
	@tex
	\beforedisplay
	$$
	x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
	$$
	\afterdisplay
	@end tex
	@ifinfo

	@example
	x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
	@end example

	@end ifinfo
	@noindent
	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, @math{y = x'(t)},
	@tex
	\beforedisplay
	$$
	\eqalign{
	x' &= y\cr
	y' &= -x + \mu y (1-x^2)
	}
	$$
	\afterdisplay
	@end tex
	@ifinfo

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

	@end ifinfo
	@noindent
	The program begins by defining functions for these derivatives and
	their Jacobian,

	@example
	@verbatiminclude examples/ode-initval.c
	@end example

	@noindent
	For functions with multiple parameters, the appropriate information
	can be passed in through the @var{params} argument using a pointer to
	a struct.

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

	@iftex
	@sp 1
	@center @image{vdp,3.4in}
	@center Numerical solution of the Van der Pol oscillator equation
	@center using Prince-Dormand 8th order Runge-Kutta.
	@end iftex

	@noindent
	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 @math{t = 0,
	1, 2, \dots, 100},

	@example
	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]);
	@}
	@end example

	@noindent
	It is also possible to work with a non-adaptive integrator, using only
	the stepping function itself. The following program uses the @code{rk4}
	fourth-order Runge-Kutta stepping function with a fixed stepsize of
	0.01,

	@example
	@verbatiminclude examples/odefixed.c
	@end example

	@noindent
	The derivatives must be initialized for the starting point @math{t=0}
	before the first step is taken. Subsequent steps use the output
	derivatives @var{dydt_out} as inputs to the next step by copying their
	values into @var{dydt_in}.

	@node ODE References and Further Reading
	@section References and Further Reading

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

	@itemize @asis
	@item
	Abramowitz & Stegun (eds.), @cite{Handbook of Mathematical Functions},
	Section 25.5.
	@end itemize

	@noindent
	The implicit Bulirsch-Stoer algorithm @code{bsimp} is described in the
	following paper,

	@itemize @asis
	@item
	G. Bader and P. Deuflhard, ``A Semi-Implicit Mid-Point Rule for Stiff
	Systems of Ordinary Differential Equations.'', Numer.@: Math.@: 41, 373--398,
	1983.
	@end itemize