.. index::

   single: differential equations, initial value problems

   single: initial value problems, differential equations

   single: ordinary differential equations, initial value problem

   single: ODEs, initial value problems

*******************************

Ordinary Differential Equations

*******************************

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. A driver object

can be used as a high level wrapper for easy use of low level

functions.

These functions are declared in the header file :file:`gsl_odeiv2.h`.

This is a new interface in version 1.15 and uses the prefix

:code:`gsl_odeiv2` for all functions.  It is recommended over the

previous :code:`gsl_odeiv` implementation defined in :file:`gsl_odeiv.h`

The old interface has been retained under the original name for

backwards compatibility.

Defining the ODE System

=======================

The routines solve the general :math:`n`-dimensional first-order system,

.. only:: not texinfo

   .. math:: {dy_i(t) \over dt} = f_i (t, y_1(t), \dots y_n(t))

.. only:: texinfo

   ::

      dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t))

for :math:`i = 1, \dots, n`.  The stepping functions rely on the vector

of derivatives :math:`f_i` and the Jacobian matrix,

.. only:: not texinfo

   .. math:: J_{ij} = \partial f_i(t, y(t)) / \partial y_j

.. only:: texinfo

   ::

      J_{ij} = df_i(t,y(t)) / dy_j

A system of equations is defined using the :type:`gsl_odeiv2_system`

datatype.

.. type:: gsl_odeiv2_system

   This data type defines a general ODE system with arbitrary parameters.

   :code:`int (* function) (double t, const double y[], double dydt[], void * params)`

      This function should store the vector elements

      :math:`f_i(t,y,params)` in the array :data:`dydt`,

      for arguments (:data:`t`, :data:`y`) and parameters :data:`params`.

      The function should return :macro:`GSL_SUCCESS` if the calculation was

      completed successfully. Any other return value indicates an error. A

      special return value :macro:`GSL_EBADFUNC` causes :code:`gsl_odeiv2`

      routines to immediately stop and return. If :code:`function` 

      is modified (for example contents of :data:`params`), the user must call an

      appropriate reset function (:func:`gsl_odeiv2_driver_reset`,

      :func:`gsl_odeiv2_evolve_reset` or :func:`gsl_odeiv2_step_reset`)

      before continuing. Use return values

      distinct from standard GSL error codes to distinguish your function as

      the source of the error.

   .. index::

      single: Jacobian matrix, ODEs

   :code:`int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params)`

      This function should store the vector of derivative elements

      .. only:: not texinfo

         .. math:: \partial f_i(t,y,params) / \partial t

      .. only:: texinfo

         ::

            df_i(t,y,params)/dt

      in the array :data:`dfdt` and the Jacobian matrix :math:`J_{ij}`

      in the array :data:`dfdy`, regarded as a row-ordered

      matrix :code:`J(i,j) = dfdy[i * dimension + j]` where :code:`dimension`

      is the dimension of the system. 

      Not all of the stepper algorithms of :code:`gsl_odeiv2` make use of the

      Jacobian matrix, so it may not be necessary to provide this function

      (the :code:`jacobian` element of the struct can be replaced by a null

      pointer for those algorithms).

      The function should return :macro:`GSL_SUCCESS` if the calculation was

      completed successfully. Any other return value indicates an error. A

      special return value :macro:`GSL_EBADFUNC` causes :code:`gsl_odeiv2`

      routines to immediately stop and return. If :code:`jacobian` 

      is modified (for example contents of :data:`params`), the user must call an

      appropriate reset function (:func:`gsl_odeiv2_driver_reset`,

      :func:`gsl_odeiv2_evolve_reset` or :func:`gsl_odeiv2_step_reset`)

      before continuing. Use return values

      distinct from standard GSL error codes to distinguish your function as

      the source of the error.

   :code:`size_t dimension`

      This is the dimension of the system of equations.

   :code:`void * params`

      This is a pointer to the arbitrary parameters of the system.

Stepping Functions

==================

The lowest level components are the *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.

.. type:: gsl_odeiv2_step

   This contains internal parameters for a stepping function.

.. function:: gsl_odeiv2_step * gsl_odeiv2_step_alloc (const gsl_odeiv2_step_type * T, size_t dim)

   This function returns a pointer to a newly allocated instance of a

   stepping function of type :data:`T` for a system of :data:`dim`

   dimensions. Please note that if you use a stepper method that

   requires access to a driver object, it is advisable to use a driver

   allocation method, which automatically allocates a stepper, too.

.. function:: int gsl_odeiv2_step_reset (gsl_odeiv2_step * s)

   This function resets the stepping function :data:`s`.  It should be used

   whenever the next use of :data:`s` will not be a continuation of a

   previous step.

.. function:: void gsl_odeiv2_step_free (gsl_odeiv2_step * s)

   This function frees all the memory associated with the stepping function

   :data:`s`.

.. function:: const char * gsl_odeiv2_step_name (const gsl_odeiv2_step * s)

   This function returns a pointer to the name of the stepping function.

   For example::

      printf ("step method is '%s'\n", gsl_odeiv2_step_name (s));

   would print something like :code:`step method is 'rkf45'`.

.. function:: unsigned int gsl_odeiv2_step_order (const gsl_odeiv2_step * s)

   This function returns the order of the stepping function on the previous

   step. The order can vary if the stepping function itself is adaptive.

.. function:: int gsl_odeiv2_step_set_driver (gsl_odeiv2_step * s, const gsl_odeiv2_driver * d)

   This function sets a pointer of the driver object :data:`d` for stepper

   :data:`s`, to allow the stepper to access control (and evolve) object

   through the driver object. This is a requirement for some steppers, to

   get the desired error level for internal iteration of

   stepper. Allocation of a driver object calls this function

   automatically.

.. function:: int gsl_odeiv2_step_apply (gsl_odeiv2_step * s, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv2_system * sys)

   This function applies the stepping function :data:`s` to the system of

   equations defined by :data:`sys`, using the step-size :data:`h` to advance

   the system from time :data:`t` and state :data:`y` to time :data:`t` + :data:`h`.

   The new state of the system is stored in :data:`y` on output, with an

   estimate of the absolute error in each component stored in :data:`yerr`.

   If the argument :data:`dydt_in` is not null it should point an array

   containing the derivatives for the system at time :data:`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 :data:`t` + :data:`h` will

   be stored in :data:`dydt_out` if it is not null.

   The stepping function returns :macro:`GSL_FAILURE` if it is unable to

   compute the requested step. Also, if the user-supplied functions

   defined in the system :data:`sys` return a status other than

   :macro:`GSL_SUCCESS` the step will be aborted. In that case, the

   elements of :data:`y` will be restored to their pre-step values and the

   error code from the user-supplied function will be returned. Failure

   may be due to a singularity in the system or too large step-size

   :data:`h`. In that case the step should be attempted again with a

   smaller step-size, e.g. :data:`h` / 2.

   If the driver object is not appropriately set via

   :func:`gsl_odeiv2_step_set_driver` for those steppers that need it, the

   stepping function returns :macro:`GSL_EFAULT`. If the user-supplied

   functions defined in the system :data:`sys` returns :macro:`GSL_EBADFUNC`,

   the function returns immediately with the same return code. In this

   case the user must call :func:`gsl_odeiv2_step_reset` before calling

   this function again.

The following algorithms are available,

.. type:: gsl_odeiv2_step_type

   .. index::

      single: RK2, Runge-Kutta method

      single: Runge-Kutta methods, ordinary differential equations

   .. var:: gsl_odeiv2_step_rk2

      Explicit embedded Runge-Kutta (2, 3) method.

   .. index::

      single: RK4, Runge-Kutta method

   .. var:: gsl_odeiv2_step_rk4

      Explicit 4th order (classical) Runge-Kutta. Error estimation is

      carried out by the step doubling method. For more efficient estimate

      of the error, use the embedded methods described below.

   .. index::

      single: Fehlberg method, differential equations

      single: RKF45, Runge-Kutta-Fehlberg method

   .. var:: gsl_odeiv2_step_rkf45

      Explicit embedded Runge-Kutta-Fehlberg (4, 5) method.  This method is

      a good general-purpose integrator.

   .. index::

      single: Runge-Kutta Cash-Karp method

      single: Cash-Karp, Runge-Kutta method

   .. var:: gsl_odeiv2_step_rkck 

      Explicit embedded Runge-Kutta Cash-Karp (4, 5) method.

   .. index::

      single: Runge-Kutta Prince-Dormand method

      single: Prince-Dormand, Runge-Kutta method

   .. var:: gsl_odeiv2_step_rk8pd  

      Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method.

   .. index:: Implicit Euler method

   .. var:: gsl_odeiv2_step_rk1imp 

      Implicit Gaussian first order Runge-Kutta. Also known as implicit

      Euler or backward Euler method. Error estimation is carried out by the

      step doubling method. This algorithm requires the Jacobian and 

      access to the driver object via :func:`gsl_odeiv2_step_set_driver`.

   .. index:: Implicit Runge-Kutta method

   .. var:: gsl_odeiv2_step_rk2imp  

      Implicit Gaussian second order Runge-Kutta. Also known as implicit

      mid-point rule. Error estimation is carried out by the step doubling

      method. This stepper requires the Jacobian and access to the driver

      object via :func:`gsl_odeiv2_step_set_driver`.

   .. var:: gsl_odeiv2_step_rk4imp  

      Implicit Gaussian 4th order Runge-Kutta. Error estimation is carried

      out by the step doubling method. This algorithm requires the Jacobian

      and access to the driver object via :func:`gsl_odeiv2_step_set_driver`.

   .. index::

      single: Bulirsch-Stoer method

      single: Bader and Deuflhard, Bulirsch-Stoer method.

      single: Deuflhard and Bader, Bulirsch-Stoer method.

   .. var:: gsl_odeiv2_step_bsimp   

      Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is

      generally suitable for stiff problems. This stepper requires the

      Jacobian.

   .. index::

      single: Adams method

      single: multistep methods, ODEs

      single: predictor-corrector method, ODEs

      single: Nordsieck form

   .. var:: gsl_odeiv2_step_msadams  

      A variable-coefficient linear multistep Adams method in Nordsieck

      form. This stepper uses explicit Adams-Bashforth (predictor) and

      implicit Adams-Moulton (corrector) methods in :math:`P(EC)^m`

      functional iteration mode. Method order varies dynamically between 1

      and 12. This stepper requires the access to the driver object via

      :func:`gsl_odeiv2_step_set_driver`.

   .. index:: BDF method

   .. var:: gsl_odeiv2_step_msbdf

      A variable-coefficient linear multistep backward differentiation

      formula (BDF) method in Nordsieck form. This stepper uses the explicit

      BDF formula as predictor and implicit BDF formula as corrector. A

      modified Newton iteration method is used to solve the system of

      non-linear equations. Method order varies dynamically between 1 and

      5. The method is generally suitable for stiff problems. This stepper

      requires the Jacobian and the access to the driver object via

      :func:`gsl_odeiv2_step_set_driver`.

.. index::

   single: Adaptive step-size control, differential equations

Adaptive Step-size Control

==========================

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.

.. type:: gsl_odeiv2_control

   This is a workspace for controlling step size.

.. type:: gsl_odeiv2_control_type

   This specifies the control type.

.. function:: gsl_odeiv2_control * gsl_odeiv2_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 :data:`eps_abs` and :data:`eps_rel`, and

   scaling factors :data:`a_y` and :data:`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,

   .. only:: not texinfo

      .. math:: D_i = \epsilon_{abs} + \epsilon_{rel} * (a_{y} |y_i| + a_{dydt} h |y\prime_i|)

   .. only:: texinfo

      ::

         D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y\prime_i|)

   and comparing it with the observed error :math:`E_i = |yerr_i|`.  If the

   observed error :data:`E` exceeds the desired error level :data:`D` by more

   than 10% for any component then the method reduces the step-size by an

   appropriate factor,

   .. math:: h_{new} = h_{old} * S * (E/D)^{-1/q}

   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 :data:`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,

   .. math:: h_{new} = h_{old} * S * (E/D)^{-1/(q+1)}

   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.

.. function:: gsl_odeiv2_control * gsl_odeiv2_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 :data:`eps_abs` and

   relative error of :data:`eps_rel` with respect to the solution :math:`y_i(t)`.

   This is equivalent to the standard control object with :data:`a_y` = 1 and

   :data:`a_dydt` = 0.

.. function:: gsl_odeiv2_control * gsl_odeiv2_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 :data:`eps_abs` and

   relative error of :data:`eps_rel` with respect to the derivatives of the

   solution :math:`y'_i(t)`.  This is equivalent to the standard control

   object with :data:`a_y` = 0 and :data:`a_dydt` = 1.

.. function:: gsl_odeiv2_control * gsl_odeiv2_control_scaled_new (double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim)

   This function creates a new control object which uses the same algorithm

   as :func:`gsl_odeiv2_control_standard_new` but with an absolute error

   which is scaled for each component by the array :data:`scale_abs`.

   The formula for :math:`D_i` for this control object is,

   .. only:: not texinfo

      .. math:: D_i = \epsilon_{abs} s_i + \epsilon_{rel} * (a_{y} |y_i| + a_{dydt} h |y\prime_i|)

   .. only:: texinfo

      ::

         D_i = eps_abs * s_i + eps_rel * (a_y |y_i| + a_dydt h |y\prime_i|)

   where :math:`s_i` is the :math:`i`-th component of the array :data:`scale_abs`.

   The same error control heuristic is used by the Matlab ODE suite. 

.. function:: gsl_odeiv2_control * gsl_odeiv2_control_alloc (const gsl_odeiv2_control_type * T)

   This function returns a pointer to a newly allocated instance of a

   control function of type :data:`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_odeiv2_control_init (gsl_odeiv2_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt)

   This function initializes the control function :data:`c` with the

   parameters :data:`eps_abs` (absolute error), :data:`eps_rel` (relative

   error), :data:`a_y` (scaling factor for y) and :data:`a_dydt` (scaling

   factor for derivatives).

.. function:: void gsl_odeiv2_control_free (gsl_odeiv2_control * c)

   This function frees all the memory associated with the control function

   :data:`c`.

.. function:: int gsl_odeiv2_control_hadjust (gsl_odeiv2_control * c, gsl_odeiv2_step * s, const double y[], const double yerr[], const double dydt[], double * h)

   This function adjusts the step-size :data:`h` using the control function

   :data:`c`, and the current values of :data:`y`, :data:`yerr` and :data:`dydt`.

   The stepping function :data:`step` is also needed to determine the order

   of the method.  If the error in the y-values :data:`yerr` is found to be

   too large then the step-size :data:`h` is reduced and the function returns

   :macro:`GSL_ODEIV_HADJ_DEC`.  If the error is sufficiently small then

   :data:`h` may be increased and :macro:`GSL_ODEIV_HADJ_INC` is returned.  The

   function returns :macro:`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_odeiv2_control_name (const gsl_odeiv2_control * c)

   This function returns a pointer to the name of the control function.

   For example::

      printf ("control method is '%s'\n", gsl_odeiv2_control_name (c));

   would print something like :code:`control method is 'standard'`

.. function:: int gsl_odeiv2_control_errlevel (gsl_odeiv2_control * c, const double y, const double dydt, const double h, const size_t ind, double * errlev)

   This function calculates the desired error level of the :data:`ind`-th component

   to :data:`errlev`. It requires the value (:data:`y`) and value of the derivative

   (:data:`dydt`) of the component, and the current step size :data:`h`.

.. function:: int gsl_odeiv2_control_set_driver (gsl_odeiv2_control * c, const gsl_odeiv2_driver * d)

   This function sets a pointer of the driver object :data:`d` for control

   object :data:`c`.

Evolution

=========

The evolution function combines the results of a stepping function and

control function to reliably advance the solution forward one step

using an acceptable step-size.

.. type:: gsl_odeiv2_evolve

   This workspace contains parameters for controlling the evolution function

.. function:: gsl_odeiv2_evolve * gsl_odeiv2_evolve_alloc (size_t dim)

   This function returns a pointer to a newly allocated instance of an

   evolution function for a system of :data:`dim` dimensions.

.. function:: int gsl_odeiv2_evolve_apply (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con, gsl_odeiv2_step * step, const gsl_odeiv2_system * sys, double * t, double t1, double * h, double y[])

   This function advances the system (:data:`e`, :data:`sys`) from time

   :data:`t` and position :data:`y` using the stepping function :data:`step`.

   The new time and position are stored in :data:`t` and :data:`y` on output.

   The initial step-size is taken as :data:`h`. The control function

   :data:`con` is applied to check whether the local error estimated by the

   stepping function :data:`step` using step-size :data:`h` exceeds the

   required error tolerance. If the error is too high, the step is

   retried by calling :data:`step` with a decreased step-size. This process

   is continued until an acceptable step-size is found. An estimate of

   the local error for the step can be obtained from the components of

   the array :code:`e->yerr[]`.

   If the user-supplied functions defined in the system :data:`sys` returns

   :macro:`GSL_EBADFUNC`, the function returns immediately with the same

   return code. In this case the user must call

   :func:`gsl_odeiv2_step_reset` and

   :func:`gsl_odeiv2_evolve_reset` before calling this function again.

   Otherwise, if the user-supplied functions defined in the system

   :data:`sys` or the stepping function :data:`step` return a status other

   than :macro:`GSL_SUCCESS`, the step is retried with a decreased

   step-size. If the step-size decreases below machine precision, a

   status of :macro:`GSL_FAILURE` is returned if the user functions

   returned :macro:`GSL_SUCCESS`. Otherwise the value returned by user

   function is returned. If no acceptable step can be made, :data:`t` and

   :data:`y` will be restored to their pre-step values and :data:`h` contains

   the final attempted step-size.

   If the step is successful the function returns a suggested step-size

   for the next step in :data:`h`. The maximum time :data:`t1` is guaranteed

   not to be exceeded by the time-step. On the final time-step the value

   of :data:`t` will be set to :data:`t1` exactly.

.. function:: int gsl_odeiv2_evolve_apply_fixed_step (gsl_odeiv2_evolve * e, gsl_odeiv2_control * con, gsl_odeiv2_step * step, const gsl_odeiv2_system * sys, double * t, const double h, double y[])

   This function advances the ODE-system (:data:`e`, :data:`sys`, :data:`con`)

   from time :data:`t` and position :data:`y` using the stepping function

   :data:`step` by a specified step size :data:`h`. If the local error

   estimated by the stepping function exceeds the desired error level,

   the step is not taken and the function returns

   :macro:`GSL_FAILURE`. Otherwise the value returned by user function is

   returned.

.. function:: int gsl_odeiv2_evolve_reset (gsl_odeiv2_evolve * e)

   This function resets the evolution function :data:`e`.  It should be used

   whenever the next use of :data:`e` will not be a continuation of a

   previous step.

.. function:: void gsl_odeiv2_evolve_free (gsl_odeiv2_evolve * e)

   This function frees all the memory associated with the evolution function

   :data:`e`.

.. function:: int gsl_odeiv2_evolve_set_driver (gsl_odeiv2_evolve * e, const gsl_odeiv2_driver * d)

   This function sets a pointer of the driver object :data:`d` for evolve

   object :data:`e`.

.. index::

   single: discontinuities, in ODE systems

If a system has discontinuous changes in the derivatives at known

points, it is advisable to evolve the system between each discontinuity

in sequence.  For example, if a step-change in an external driving

force occurs at times :math:`t_a, t_b` and :math:`t_c` then evolution

should be carried out over the ranges :math:`(t_0,t_a)`,

:math:`(t_a,t_b)`, :math:`(t_b,t_c)`, and :math:`(t_c,t_1)` separately

and not directly over the range :math:`(t_0,t_1)`.

Driver

======

The driver object is a high level wrapper that combines the evolution,

control and stepper objects for easy use.

.. function:: gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel)

              gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel)

              gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel, const double a_y, const double a_dydt)

              gsl_odeiv2_driver * gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys, const gsl_odeiv2_step_type * T, const double hstart, const double epsabs, const double epsrel, const double a_y, const double a_dydt, const double scale_abs[])

   These functions return a pointer to a newly allocated instance of a

   driver object. The functions automatically allocate and initialise the

   evolve, control and stepper objects for ODE system :data:`sys` using

   stepper type :data:`T`. The initial step size is given in

   :data:`hstart`. The rest of the arguments follow the syntax and

   semantics of the control functions with same name

   (:code:`gsl_odeiv2_control_*_new`).

.. function:: int gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin)

   The function sets a minimum for allowed step size :data:`hmin` for

   driver :data:`d`. Default value is 0.

.. function:: int gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax)

   The function sets a maximum for allowed step size :data:`hmax` for

   driver :data:`d`. Default value is :macro:`GSL_DBL_MAX`.

.. function:: int gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d, const unsigned long int nmax)

   The function sets a maximum for allowed number of steps :data:`nmax` for

   driver :data:`d`. Default value of 0 sets no limit for steps.

.. function:: int gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double * t, const double t1, double y[])

   This function evolves the driver system :data:`d` from :data:`t` to

   :data:`t1`. Initially vector :data:`y` should contain the values of

   dependent variables at point :data:`t`. If the function is unable to

   complete the calculation, an error code from

   :func:`gsl_odeiv2_evolve_apply` is returned, and :data:`t` and :data:`y`

   contain the values from last successful step. 

   If maximum number of steps is reached, a value of :macro:`GSL_EMAXITER`

   is returned. If the step size drops below minimum value, the function

   returns with :macro:`GSL_ENOPROG`. If the user-supplied functions

   defined in the system :data:`sys` returns :macro:`GSL_EBADFUNC`, the

   function returns immediately with the same return code. In this case

   the user must call :func:`gsl_odeiv2_driver_reset` before calling this

   function again.

.. function:: int gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double * t, const double h, const unsigned long int n, double y[])

   This function evolves the driver system :data:`d` from :data:`t` with

   :data:`n` steps of size :data:`h`. If the function is unable to complete

   the calculation, an error code from

   :func:`gsl_odeiv2_evolve_apply_fixed_step` is returned, and :data:`t` and

   :data:`y` contain the values from last successful step.

.. function:: int gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d)

   This function resets the evolution and stepper objects.

.. function:: int gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart)

   The routine resets the evolution and stepper objects and sets new

   initial step size to :data:`hstart`. This function can be used e.g. to

   change the direction of integration.

.. function:: int gsl_odeiv2_driver_free (gsl_odeiv2_driver * d)

   This function frees the driver object, and the related evolution,

   stepper and control objects.

Examples

========

.. index::

   single: Van der Pol oscillator, example

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

oscillator equation,

.. math:: u''(t) + \mu u'(t) (u(t)^2 - 1) + u(t) = 0

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:`v = u'(t)`,

.. only:: not texinfo

   .. math::

      u' &= v \\

      v' &= -u + \mu v (1-u^2)

.. only:: texinfo

   ::

      u' = v

      v' = -u + \mu v (1-u^2)

The program begins by defining functions for these derivatives and

their Jacobian. The main function uses driver level functions to solve

the problem. The program evolves the solution from :math:`(u, v) = (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 :math:`10^{-6}`

in the function values :math:`(u, v)`.

The loop in the example prints the solution at the points

:math:`t_i = 1, 2, \dots, 100`.

.. include:: examples/ode-initval.c

   :code:

The user can work with the lower level functions directly, as in

the following example. In this case an intermediate result is printed

after each successful step instead of equidistant time points. 

.. include:: examples/ode-initval-low-level.c

   :code:

For functions with multiple parameters, the appropriate information

can be passed in through the :data:`params` argument in

:type:`gsl_odeiv2_system` definition (:data:`mu` in this example) by using

a pointer to a struct.

.. figure:: /images/ode-vdp.png

   Numerical solution of the Van der Pol oscillator equation 

   using Prince-Dormand 8th order Runge-Kutta.

It is also possible to work with a non-adaptive integrator, using only

the stepping function itself,

:func:`gsl_odeiv2_driver_apply_fixed_step` or

:func:`gsl_odeiv2_evolve_apply_fixed_step`. The following program uses

the driver level function, with fourth-order

Runge-Kutta stepping function with a fixed stepsize of

0.001.

.. include:: examples/odefixed.c

   :code:

References and Further Reading

==============================

* Ascher, U.M., Petzold, L.R., *Computer Methods for Ordinary

  Differential and Differential-Algebraic Equations*, SIAM, 

  Philadelphia, 1998.

* Hairer, E., Norsett, S. P., Wanner, G., *Solving Ordinary Differential 

  Equations I: Nonstiff Problems*, Springer, Berlin, 1993.

* Hairer, E., Wanner, G., *Solving Ordinary Differential 

  Equations II: Stiff and Differential-Algebraic Problems*,

  Springer, Berlin, 1996.

Many of 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 :code:`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.

The Adams and BDF multistep methods :code:`msadams` and :code:`msbdf`

are based on the following articles,

* G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the

  Numerical Solution of Ordinary Differential Equations.",

  ACM Trans. Math. Software, 1, 71--96, 1975.

* P. N. Brown, G. D. Byrne and A. C. Hindmarsh, "VODE: A

  Variable-coefficient ODE Solver.", SIAM J. Sci. Stat. Comput. 10,

  1038--1051, 1989.

* A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban,

  D. E. Shumaker and C. S. Woodward, "SUNDIALS: Suite of

  Nonlinear and Differential/Algebraic Equation Solvers.", ACM

  Trans. Math. Software 31, 363--396, 2005.