diff options
Diffstat (limited to 'gsl-1.9/doc/ode-initval.texi')
| -rw-r--r-- | gsl-1.9/doc/ode-initval.texi | 556 |
1 files changed, 556 insertions, 0 deletions
diff --git a/gsl-1.9/doc/ode-initval.texi b/gsl-1.9/doc/ode-initval.texi new file mode 100644 index 0000000..ae8e68b --- /dev/null +++ b/gsl-1.9/doc/ode-initval.texi @@ -0,0 +1,556 @@ +@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 |