10.2 Ordinary differential equation (ODE) solvers
Stan provides several higher order functions for solving initial value problems specified as Ordinary Differential Equations (ODEs).
Solving an initial value ODE means given a set of differential equations \(y'(t, \theta) = f(t, y, \theta)\) and initial conditions \(y(t_0, \theta)\), solving for \(y\) at a sequence of times \(t_0 < t_1 < t_2, \cdots < t_n\). \(f(t, y, \theta)\) is referred to here as the ODE system function.
\(f(t, y, \theta)\) will be defined as a function with a certain signature and provided along with the initial conditions and output times to one of the ODE solver functions.
To make it easier to write ODEs, the solve functions take extra arguments
that are passed along unmodified to the user-supplied system function.
Because there can be any number of these arguments and they can be of different types,
they are denoted below as .... The types of the arguments represented by ...
in the ODE solve function call must match the types of the arguments represented by
... in the user-supplied system function.
10.2.1 Non-stiff solver
array[] vector ode_rk45(function ode, vector initial_state, real initial_time, array[] real times, ...)
Solves the ODE system for the times provided using the Dormand-Prince
algorithm, a 4th/5th order Runge-Kutta method.
Available since 2.24
array[] vector ode_rk45_tol(function ode, vector initial_state, real initial_time, array[] real times, data real rel_tol, data real abs_tol, int max_num_steps, ...)
Solves the ODE system for the times provided using the Dormand-Prince
algorithm, a 4th/5th order Runge-Kutta method with additional control
parameters for the solver.
Available since 2.24
array[] vector ode_ckrk(function ode, vector initial_state, real initial_time, array[] real times, ...)
Solves the ODE system for the times provided using the Cash-Karp
algorithm, a 4th/5th order explicit Runge-Kutta method.
Available since 2.27
array[] vector ode_ckrk_tol(function ode, vector initial_state, real initial_time, array[] real times, data real rel_tol, data real abs_tol, int max_num_steps, ...)
Solves the ODE system for the times provided using the Cash-Karp
algorithm, a 4th/5th order explicit Runge-Kutta method with additional control
parameters for the solver.
Available since 2.27
array[] vector ode_adams(function ode, vector initial_state, real initial_time, array[] real times, ...)
Solves the ODE system for the times provided using the Adams-Moulton method.
Available since 2.24
array[] vector ode_adams_tol(function ode, vector initial_state, real initial_time, array[] real times, data real rel_tol, data real abs_tol, int max_num_steps, ...)
Solves the ODE system for the times provided using the Adams-Moulton
method with additional control parameters for the solver.
Available since 2.24
10.2.2 Stiff solver
array[] vector ode_bdf(function ode, vector initial_state, real initial_time, array[] real times, ...)
Solves the ODE system for the times provided using the backward differentiation
formula (BDF) method.
Available since 2.24
array[] vector ode_bdf_tol(function ode, vector initial_state, real initial_time, array[] real times, data real rel_tol, data real abs_tol, int max_num_steps, ...)
Solves the ODE system for the times provided using the backward differentiation
formula (BDF) method with additional control parameters for the solver.
Available since 2.24
10.2.3 Adjoint solver
array[] vector ode_adjoint_tol_ctl(function ode, vector initial_state, real initial_time, array[] real times, data real rel_tol_forward, data vector abs_tol_forward, data real rel_tol_backward, data vector abs_tol_backward, int max_num_steps, int num_steps_between_checkpoints, int interpolation_polynomial, int solver_forward, int solver_backward, ...)
Solves the ODE system for the times provided using the adjoint ODE solver method
from CVODES. The adjoint ODE solver requires a checkpointed forward in time ODE
integration, a backwards in time integration that makes uses of an interpolated
version of the forward solution, and the solution of a quadrature problem (the
number of which depends on the number of parameters passed to the solve). The
tolerances and numeric methods used for the forward solve, backward solve,
quadratures, and interpolation can all be configured.
Available since 2.27
10.2.4 ODE system function
The first argument to one of the ODE solvers is always the ODE system
function. The ODE system function must have a vector return type, and the
first two arguments must be a real and vector in that order. These two
arguments are followed by the variadic arguments that are passed through from
the ODE solve function call:
vector ode(real time, vector state, ...)The ODE system function should return the derivative of the state with respect to time at the time and state provided. The length of the returned vector must match the length of the state input into the function.
The arguments to this function are:
time, the time to evaluate the ODE systemstate, the state of the ODE system at the time specified..., sequence of arguments passed unmodified from the ODE solve function call. The types here must match the types in the...arguments of the ODE solve function call.
10.2.5 Arguments to the ODE solvers
The arguments to the ODE solvers in both the stiff and non-stiff solvers are the same. The arguments to the adjoint ODE solver are different; see Arguments to the adjoint ODE solvers.
ode: ODE system function,initial_state: initial state, typevector,initial_time: initial time, typereal,times: solution times, typearray[] real,...: sequence of arguments that will be passed through unmodified to the ODE system function. The types here must match the types in the...arguments of the ODE system function.
For the versions of the ode solver functions ending in _tol, these three
parameters must be provided after times and before the ... arguments:
datarel_tol: relative tolerance for the ODE solver, typereal, data only,dataabs_tol: absolute tolerance for the ODE solver, typereal, data only, andmax_num_steps: maximum number of steps to take between output times in the ODE solver, typeint, data only.
Because the tolerances are data arguments, they must be defined in either the data
or transformed data blocks. They cannot be parameters, transformed parameters
or functions of parameters or transformed parameters.
10.2.6 Arguments to the adjoint ODE solver
The arguments to the adjoint ODE solver are different from those for the other functions (for those see Arguments to the adjoint ODE solvers).
ode: ODE system function,initial_state: initial state, typevector,initial_time: initial time, typereal,times: solution times, typearray[] real,datarel_tol_forward: Relative tolerance for forward solve, typereal, data only,dataabs_tol_forward: Absolute tolerance vector for each state for forward solve, typevector, data only,datarel_tol_backward: Relative tolerance for backward solve, typereal, data only,dataabs_tol_backward: Absolute tolerance vector for each state for backward solve, typevector, data only,datarel_tol_quadrature: Relative tolerance for backward quadrature, typereal, data only,dataabs_tol_quadrature: Absolute tolerance for backward quadrature, typereal, data only,datamax_num_steps: Maximum number of time-steps to take in integrating the ODE solution between output time points for forward and backward solve, typeint, data only,num_steps_between_checkpoints: number of steps between checkpointing forward solution, typeint, data only,interpolation_polynomial: can be 1 for hermite or 2 for polynomial interpolation method of CVODES, typeint, data only,solver_forward: solver used for forward ODE problem: 1=Adams (non-stiff), 2=BDF (stiff), typeint, data only,solver_backward: solver used for backward ODE problem: 1=Adams (non-stiff), 2=BDF (stiff), typeint, data only....: sequence of arguments that will be passed through unmodified to the ODE system function. The types here must match the types in the...arguments of the ODE system function.
Because the tolerances are data arguments, they must be defined in either the data
or transformed data blocks. They cannot be parameters, transformed parameters
or functions of parameters or transformed parameters.
10.2.6.1 Return values
The return value for the ODE solvers is an array of vectors (type array[] vector),
one vector representing the state of the system at every time in specified in
the times argument.
10.2.6.2 Array and vector sizes
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 length of each vector in the output,
number of solution times and number of vectors in the output.