This is an old version, view current version.

13.2 Coding an ODE System

A system of ODEs is coded directly in Stan as a function with a strictly specified signature. For example, the simple harmonic oscillator can be coded using the following function in Stan (see the user-defined functions chapter for more information on coding user-defined functions).

real[] sho(real t,        // time
           real[] y,      // state
           real[] theta,  // parameters
           real[] x_r,    // data (real)
           int[] x_i) {   // data (integer)
  real dydt[2];
  dydt[1] = y[2];
  dydt[2] = -y[1] - theta[1] * y[2];
  return dydt;
}

The function takes in a time t (a real value), a a system state y (a real array), system parameters theta (a real array), along with real data in variable x_r (a real array) and integer data in variable x_i (an integer array). The system function returns the array of derivatives of the system state with respect to time, evaluated at time t and state y. The simple harmonic oscillator coded here does not have time-sensitive equations; that is, t does not show up in the definition of dydt. The simple harmonic oscillator does not use real or integer data, either. Nevertheless, these unused arguments must be included as arguments in the system function with exactly the signature shown above.

Strict Signature

The function defining the system must have exactly these argument types and return type. This may require passing in zero-length arrays for data or parameters if the system does not involve data or parameters. A full example for the simple harmonic oscillator, which does not depend on any constant data variables, is provided in the simple harmonic oscillator trajectory plot.

Discontinuous ODE System Function

The ODE integrator is able to integrate over discontinuities in the state function, although the accuracy of points near the discontinuity may be problematic (requiring many small steps). An example of such a discontinuity is a lag in a pharmacokinetic model, where a concentration is going to be zero for times \(0 < t < t'\) for some lag-time \(t'\), whereas it will be nonzero for times \(t \geq t'\). In this example example, we would use code in the system such as

if (t < t_lag)
  return 0;
else
  ... return non-zero value...;