.. highlight:: lua

.. module:: ode

Ordinary Differential Equations
===============================

Overview
--------

This chapter describes functions for solving ordinary differential equation (ODE) initial value problems.

.. note::
  The ODE integration routines in LGSL are based on GSL routines, but they are completely rewritten in Lua.
  Only a few integration methods are actually available for the moment such as Runge-Kutta-Fehlberg and Price-Dormand methods.

  The actual interface is also likely to be changed in the near future to handle ODE systems in array form.

In LGSL, an ODE system is integrated by using an ODE solver object.
The ODE solver internally stores the state of the solver, and you can advance the solution step-by-step until, eventually, the desired value of `t` is reached.

Creating a ODE system solver
----------------------------

An ODE solver allows users to obtain a numerical solution of an Ordinary
Differential Equation (ODE) system. The ODE solver lets you solve the
general n-dimensional first-order system,

.. math::
     \frac{\textrm{d}y_i(t)}{\textrm{d}t} = 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, for some methods, also on the Jacobian matrix,

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

For the moment, none of the methods implemented in LGSL use the Jacobian matrix.

.. note::
   The current implementation is limited to systems with a few number of variables.
   It is not recommended to use the current implementation if you have more than 20 variables.
   An implementation for ODE systems in array form should be available in the near future.

ODE solver usage example
------------------------

The following example illustrates the usage of an ODE solver. The differential equation that we want to integrate is:

   .. math::
      \begin{array}{ll}
         x' = & -y - x^2 \\
	 y' = & 2 x - y^3
      \end{array}

and this is the code that we can write to implement it::

   ode = require("lgsl.ode")
   -- define the ODE function
   function odef(t, x, y)
      return -y-x^2, 2*x - y^3
   end

   -- create the ODE solver
   s = ode{N= 2, eps_abs= 1e-8}

   -- we define initial values
   t0, t1, h0 = 0, 30, 0.04
   x0, y0 = 1, 1

   -- we initialize the ODE solver
   s:init(t0, h0, odef, x0, y0)

   -- the ODE solver is iterated till the time t1 is reached
   while s.t < t1 do
      s:step(t1)
   end

Alternatively, you may want to make a plot of the curve that you obtain. This requires the package :ref:`graph-toolkit <graphics-chapter>`. In the following example, we create a "path" to describe the curve that we want to plot. We iterate with the ODE solver and we add all points with the "line_to" method. Then we create an empty plot and we add the line that we have just created::

   graph = require("graph")

   -- we create a line and add the points obtained by integrating the ODE
   ln = graph.path(x0, y0)
   while s.t < t1 do
      s:step(t1)
      ln:line_to(s.y[1], s.y[2])
   end

   -- we create the plot by adding the line
   p = graph.plot('ODE integration example')
   p:addline(ln)
   p:show()

This the plot with the curve obtained by integration of the above ODE system:

.. figure:: ode-integration-quasi-spiral.png


A Slightly Improved Example
---------------------------

In the example given above we have used the :meth:`~ODE.step` method to advance the ODE integrator.
When you use :meth:`ODE.step`, the ODE integrator will adapt at each step the step size in order to respect the maximum absolute and relative error that you requested.
This is a quite convenient behaviour but it can have a drawback since the sampling points can be very tightly packed or very largely spaced depending on the ODE system and the integration method.
In some cases this is undesirable and you may want to obtain the values with a fixed sampling size.
In this case you can use the :meth:`~ODE.evolve` method to obtain the values with a given sampling step.

So, in the example above you can change the while loop and use the :meth:`~ODE.evolve` method instead::

   -- we create a line and add the points obtained by integrating the ODE
   ln = graph.path(x0, y0)
   for t, y1, y2 in s:evolve(t1, 0.1) do
      ln:line_to(y1, y2)
   end

to obtain values sampled with a spacing of 0.1 for the t values.

You can see that the :meth:`~ODE.evolve` method works actually as a Lua iterator.

You may also note that the :meth:`~ODE.evolve` iterator provides at each iteration the value t and each of the system variables in the standard order.


ODE Solver Class Definition
---------------------------

.. class:: ODE

   Solver for an ODE system.

   .. function:: ode(spec)

      Create a new solver for an ODE system. The ``spec`` should be a table
      containing the following fields:

      N
          The dimension or number of variables of the ODE system.
      eps_abs
          The maximum absolute error in the solution ``y`` that should be tolerated.
      eps_rel, *optional*
          The maximum relative error in the solution ``y`` that should be tolerated.
      method, *optional*
          The low-level integration method used. Can be chosen between:

          - ``"rkf45"``, Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.

          - ``"rk8pd"``, Embedded Runge-Kutta Prince-Dormand (8,9) method.

   .. method:: init(t0, h0, f, y0_1, y0_2, ..., y0_N)

      Initialize the state of the solver to the time ``t0``.
      The second argument ``h0`` is the initial step size that the integrator will try.
      The function ``f`` is the function that defines the ODE system.
      It will be called as follows: ``f(t, y_1, y_2, ..., y_N)`` where ``t`` is the time and ``y_1, y_2, ...`` are the values of the N independent variables conventionally denoted here by 'y'.
      The function ``f`` should return N values that correspond to values ``f_i(t, y_1, ..., y_N)`` for each component ``f_i`` of the ODE system function.

   .. method:: step(t1)

      Advance the solution of the system by a step chosen adaptively based on the previous step size.
      The new values (t, y) are stored internally by the solver and can be retrieved as properties with the name ``t`` and ``y`` where the latter is a column matrix of size N.
      The new values of ``t`` will be less than or equal to the value given ``t1``.
      If the value ``s.t`` is less then ``t1``, then the function evolve should be called again by the user.

   .. method:: evolve(t1, t_step)

      Returns a Lua iterator that advances the ODE system at each iteration of a step ``t_step`` until the value ``t1`` is reached.
      The iterators returns the value ``t`` itself and all the system variables ``y0``, ``y1``, ... in the standard order.

      Example::

         -- we suppose an ODE system 's' is already defined and initialized
         t1 = 50
         for t, y1, y2 in s:evolve(t1, 0.5) do
            print(t, y1, y2)
         end