LibRK4

A simple step of Runge-Kutta 4 integration.

Description

The library implements a Runge-Kutta 4 scheme with the following tableau:

      0  |  0   0   0   0
     1/2 | 1/2  0   0   0
     1/2 |  0  1/2  0   0
      1  |  0   0   1   0
     ----+----------------
         | 1/6 1/3 1/3 1/6

For the ode in the form:

The integration scheme evaluates:

where:

The input u(t) is assumed constant between integration ministeps. From an implementation point of view, p is a vector of vectors. This is done for compatibility with MATLAB System Identification Toolbox model file. Since we usually don't have an initialization and a termination callback, everything must be handled inside the integrator step (malloc and free included).

Examples

Ode 1st order:

To represent an ODE in the form:

it is necessary to define the following vector field:

void ode(double *xdot,     /**< Output of the vector field */
         double t,         /**< Current time */
         const double *x,  /**< Current state vector */
         const double *u,  /**< Input vecotr */
         const double **p, /**< Vector of vectors of parameters */
         void *data)       /**< Space for user data */
{
  xdot[0] = -p[0][0] * x[0] + p[1][0] * u[0];
}

and the options for the integrator:

rk4_opts options = {
  1e-3, /**< Time steps */
  1,    /**< vector field dimensions */
  ode   /**< vector field callback (to evaluate it) */
};

a single integration step can be performed as follows:

rk4(&options, xp, t, x, u, p, NULL);

The last NULL is userdata placeholder (the callback does not use it). The complete code is in test_ode_1. The result of the integration:

Ode 2nd order:

To represent an ODE in the form:

it is necessary to define the following vector field:

void ode(double *xdot,     /**< Output of the vector field */
         double t,         /**< Current time */
         const double *x,  /**< Current state vector */
         const double *u,  /**< Input vecotr */
         const double **p, /**< Vector of vectors of parameters */
         void *data)       /**< Space for user data */
{
  xdot[0] = x[1];
  xdot[1] = -p[0][0] * x[0] - p[0][1] * x[1] + p[1][0] * u[0];
}

and the options for the integrator:

rk4_opts options = {
  1e-3, /**< Time steps */
  2,    /**< vector field dimensions */
  ode   /**< vector field callback (to evaluate it) */
};

a single integration step can be performed as follows:

rk4(&options, xp, t, x, u, p, NULL);

The complete code is in test_ode_2. The result of the integration:

Non-linear ODE

To represent an ODE in the form:

it is necessary to define the following vector field (in this case we are evaluating the same ODE with 3 different initial conditions):

void ode(double *xdot,     /**< Output of the vector field */
         double t,         /**< Current time */
         const double *x,  /**< Current state vector */
         const double *u,  /**< Input vecotr */
         const double **p, /**< Vector of vectors of parameters */
         void *data)       /**< Space for user data */
{
  for (size_t i = 0; i < 3; i++)
    xdot[i] = (1 - 2 * t) * x[i] * x[i];
}

and the options for the integrator:

rk4_opts options = {
  1e-5, /**< Time steps */
  3,    /**< vector field dimensions */
  ode   /**< vector field callback (to evaluate it) */
};

a single integration step can be performed as follows (there is no input and no parameter):

rk4(&options, xp, t, x, NULL, NULL, NULL);

The complete code is in test_ode_3. The result of the integration: