This is an old version, view current version.

9.2 Ordinary Differential Equation (ODE) Solvers

Stan provides built-in ordinary differential equation (ODE) solvers. Although they look like function applications, the ODE solvers are special in two ways.

First, the first argument to each of the solvers is a function specifying the ODE system as an argument, like PKBugs (Lunn et al. 1999). Ordinary Stan functions do not allow functions as arguments.

Second, some of the arguments to the ODE solvers are restricted to data only expressions. These expressions must not contain variables other than those declared in the data or transformed data blocks. Ordinary Stan functions place no restriction on the origin of variables in their argument expressions.

9.2.1 Specifying an Ordinary Differential Equation as a Function

A system of ODEs is specified as an ordinary function in Stan within the functions block. The ODE system function must have this function signature:

 real[] ode(real time, real[] state, real[] theta,
            real[] x_r, int[] x_i)

The ODE system function should return the derivative of the state with respect to time at the time provided. The length of the returned real array must match the length of the state input into the function.

The arguments to this function are:

  • time, the time to evaluate the ODE system

  • state, the state of the ODE system at the time specified

  • theta, parameter values used to evaluate the ODE system

  • x_r, data values used to evaluate the ODE system

  • x_i, integer data values used to evaluate the ODE system.

The ODE system function separates parameter values, theta, from data values, x_r, for efficiency in computing the gradients of the ODE.

9.2.2 Non-Stiff Solver

real[ , ] integrate_ode_rk45(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, real[] x_r, int[] x_i)
Solves the ODE system for the times provided using the Runge Kutta Dopri algorithm with the implementation from Boost.

real[ , ] integrate_ode_rk45(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, real[] x_r, int[] x_i, real rel_tol, real abs_tol, int max_num_steps)
Solves the ODE system for the times provided using the Runge Kutta Dopri algorithm with the implementation from Boost with additional control parameters for the solver.

real[ , ] integrate_ode(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, real[] x_r, int[] x_i)
Deprecated. Solves the ODE system for the times provided with a non-stiff solver. This calls the Runge Kutta Dopri algorithm.

real[] integrate_ode_adams(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, data real[] x_r, data int[] x_i)
Solves the ODE system for the times provided using the Adams-Moulton method with the implementation from CVODES.

real[] integrate_ode_adams(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, data real[] x_r, data int[] x_i, data real rel_tol, data real abs_tol, dta int max_num_steps)
Solves the ODE system for the times provided using the Adams-Moulton method with the implementation from CVODES with additional control parameters for the CVODES solver.

9.2.3 Stiff Solver

real[] integrate_ode_bdf(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, data real[] x_r, data int[] x_i)
Solves the ODE system for the times provided using the backward differentiation formula (BDF) method with the implementation from CVODES.

real[] integrate_ode_bdf(function ode, real[] initial_state, real initial_time, real[] times, real[] theta, data real[] x_r, data int[] x_i, data real rel_tol, data real abs_tol, dta int max_num_steps)
Solves the ODE system for the times provided using the backward differentiation formula (BDF) method with the implementation from CVODES with additional control parameters for the CVODES solver.

9.2.4 Arguments to the ODE solvers

The arguments to the ODE solvers in both the stiff and non-stiff cases are as follows.

  • ode: function literal referring to a function specifying the system of differential equations with signature described in ode functions:
 (real, real[], real[], data real[], data int[]):real[]

The arguments represent (1) time, (2) system state, (3) parameters, (4) real data, and (5) integer data, and the return value contains the derivatives with respect to time of the state,

  • initial_state: initial state, type real[],

  • initial_time: initial time, type int or real,

  • times: solution times, type real[],

  • theta: parameters, type real[],

  • data x_r: real data, type real[], data only, and

  • data x_i: integer data, type int[], data only.

For more fine-grained control of the ODE solvers, these parameters can also be provided:

  • data rel_tol: relative tolerance for the ODE solver, type real, data only,

  • data abs_tol: absolute tolerance for the ODE solver, type real, data only, and

  • data max_num_steps: maximum number of steps to take in the ODE solver, type int, data only.

9.2.4.1 Return values

The return value for the ODE solvers is an array of type real[,], with values consisting of solutions at the specified times.

9.2.4.2 Sizes and parallel arrays

The sizes must match, and in particular, the following groups are of the same size:

  • state variables passed into the system function, derivatives returned by the system function, initial state passed into the solver, and rows of the return value of the solver,

  • solution times and number of rows of the return value of the solver,

  • parameters, real data and integer data passed to the solver will be passed to the system function

References

Lunn, D. J., J. Wakefield, A. Thomas, N. Best, and D. Spiegelhalter. 1999. PKBugs User Guide.