Cmdstan 2.24 introduces a new ODE interface intended to make it easier to specify the ODE system function by avoiding packing and unpacking schemes required with the old interface.

Stan solves for \(y(t, \theta)\) at a sequence of times \(t_1, t_2, \cdots, t_N\) in the ODE initial value problem defined by:

\[ y(t, \theta)' = f(t, y, \theta)\\ y(t = t_0, \theta) = y_0 \]

For notation, \(y(t, \theta)\) is the state, \(f(t, y, \theta)\) is the ODE system function, and \(y_0\) and \(t_0\) are the initial conditions.

The solution, \(y(t, \theta)\), is written explicitly in terms of \(\theta\) as a reminder that the solution of an ODE initial value problem can be a function of model parameters or data. This is the usual use case: an ODE is parameterized by a set of parameters that are to be estimated.

Specifying an ODE initial value problem in Stan involves writing a Stan function for the ODE system function, \(f(t, y, \theta)\). The big difference in the old and new ODE interfaces is that previously any arguments meant for the ODE system function had to be manually packed and unpacked from special argument arrays that were passed from the ODE solve function call to the ODE system function – in the new interface none of this packing and unpacking is necessary.

The new interface also uses `vector`

variables for the
state rather than `array[] real`

variables.

As an example, in the old solver interface the system function for the ODE \(y' = -\alpha y\) would be written:

```
functions {
array[] real rhs(real t, array[] real y, array[] real theta,
array[] real x_r, array[] int x_i) {
real alpha = theta[1];
array[1] real yp = {-alpha * y[1]};
return yp;
}
}
```

In the new interface the system function can be written:

```
functions {
vector rhs(real t, vector y, real alpha) {
vector[1] yp = -alpha * y;
return yp;
}
}
```

The new interface avoids any unused arguments (such as
`x_r`

, and `x_i`

in this example), and the
parameter `alpha`

can be passed directly instead of being
packed into `theta`

.

For a simple function, this does not look like much, but for more complicated models with numerous arguments of different types, the packing and unpacking is tedious and error prone. This leads to models that are difficult to debug and difficult to iterate on.

The new interface introduces six new functions:

`ode_bdf`

, `ode_adams`

, `ode_rk45`

and `ode_bdf_tol`

,`ode_adams_tol`

,
`ode_rk45_tol`

The solvers in the first columns have default tolerance settings. The solvers in the second column accept arguments for relative tolerance, absolute tolerance, and the maximum number of steps to take between output times.

This is different from the old interface where tolerances are presented through using the same function name with a few more arguments.

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.

The new `ode_bdf`

solver interface is (the interfaces for
`ode_adams`

and `ode_rk45`

are the same):

`vector[] ode_bdf(F f, vector y0, real t0, array[] real times, ...)`

The arguments are:

`f`

- ODE system function`y0`

- Initial state of the ODE`t0`

- Initial time of the ODE`times`

- Sorted array of times to which the ode will be solved (each element must be greater than t0, but times do not need to be strictly increasing)`...`

- Sequence of arguments passed unmodified to the ODE system function. There can be any number of`...`

arguments, and the`...`

arguments can be any type, but they must match the types of the corresponding`...`

arguments of`f`

.

The ODE system function should take the form:

`vector f(real t, vector y, ...)`

The arguments are:

`t`

- Time at which to evaluate the ODE system function`y`

- State at which to evaluate the ODE system function`...`

- Sequence of arguments passed unmodified from the ODE solver function call. The`...`

must match the types of the corresponding`...`

arguments of the ODE solver function call.

A call to `ode_bdf`

returns the solution of the ODE
specified by the system function (`f`

) and the initial
conditions (`y0`

and `t0`

) at the time points
given by the `times`

argument. The solution is given by an
array of vectors.

The `ode_bdf_tol`

interface is (the interfaces for
`ode_rk45_tol`

and `ode_adams_tol`

are the
same):

```
array[] vector ode_bdf_tol(F f, vector y0, real t0, array[] real times,
real rel_tol, real abs_tol, int max_num_steps, ...)
```

The arguments are: 1. `f`

- ODE system function

`y0`

- Initial state of the ODE`t0`

- Initial time of the ODE`times`

- Sorted array of times to which the ode will be solved (each element must be greater than t0, but times do not need to be strictly increasing)`rel_tol`

- Relative tolerance for solve (data)`abs_tol`

- Absolute tolerance for solve (data)`max_num_steps`

- Maximum number of timesteps to take in integrating the ODE solution between output time points (data)`...`

- Sequence of arguments passed unmodified to the ODE system function. There can be any number of`...`

arguments, and the`...`

arguments can be any type, but they must match the types of the corresponding`...`

arguments of`f`

.

The `ode_rk45`

/`ode_bdf`

/`ode_adams`

interfaces are just wrappers around the
`ode_rk45_tol`

/`ode_bdf_tol`

/`ode_adams_tol`

interfaces with defaults for `rel_tol`

, `abs_tol`

,
and `max_num_steps`

. For the RK45 solver the defaults are
\(10^{-6}\) for `rel_tol`

and `abs_tol`

and \(10^6\)
for `max_num_steps`

. For the BDF/Adams solvers the defaults
are \(10^{-10}\) for
`rel_tol`

and `abs_tol`

and \(10^8\) for `max_num_steps`

.

For more detailed information about either interface, look at the function reference guide: New interface, Old interface

The two models here come from the Stan Statistical Computation Benchmarks.

In the old SIR system function, `beta`

,
`kappa`

, `gamma`

, `xi`

, and
`delta`

, are packed into the `array[] real theta`

argument. `kappa`

isn’t actually a model parameter so it is
not clear why it is packed in with the other parameters, but it is.
Promoting `kappa`

to a parameter causes there to be more
states in the extended ODE sensitivity system (used to get gradients of
the ODE with respect to inputs). Adding states to the sensitivity system
makes the ODE harder to solve and should always be avoided. The ODE
system function looks like:

```
functions {
// theta[1] = beta, water contact rate
// theta[2] = kappa, C_{50}
// theta[3] = gamma, recovery rate
// theta[4] = xi, bacteria production rate
// theta[5] = delta, bacteria removal rate
array[] real simple_SIR(real t, array[] real y, array[] real theta,
array[] real x_r, array[] int x_i) {
array[4] real dydt;
dydt[1] = -theta[1] * y[4] / (y[4] + theta[2]) * y[1];
dydt[2] = theta[1] * y[4] / (y[4] + theta[2]) * y[1] - theta[3] * y[2];
dydt[3] = theta[3] * y[2];
dydt[4] = theta[4] * y[2] - theta[5] * y[4];
return dydt;
}
}
```

For comparison, with the new interface the ODE system function can be
rewritten to explicitly name all the parameters. No separation of data
and parameters is necessary either – the solver will not add more
equations for arguments that are defined in the `data`

and
`transformed data`

blocks. The state variables in the new
model are also represented by `vector`

variables instead of
`array[]`

variables. The new ODE system function is:

```
functions {
vector simple_SIR(real t,
vector y,
real beta, // water contact rate
real kappa, // C_{50}
real gamma, // recovery rate
real xi, // bacteria production rate
real delta) { // bacteria removal rate
vector[4] dydt;
dydt[1] = -beta * y[4] / (y[4] + kappa) * y[1];
dydt[2] = beta * y[4] / (y[4] + kappa) * y[1] - gamma * y[2];
dydt[3] = gamma * y[2];
dydt[4] = xi * y[2] - delta * y[4];
return dydt;
}
}
```

In the old ODE interface, the parameters are all packed into a
`array[] real`

before calling the ODE solver:

```
transformed parameters {
array[N_t, 4] real<lower=0> y;
{
array[5] real theta = {beta, kappa, gamma, xi, delta};
y = integrate_ode_rk45(simple_SIR, y0, t0, t, theta, x_r, x_i);
}
}
```

In the new ODE interface each of the arguments is appended on to the
ODE solver function call. The RK45 ODE solver with default tolerances is
called `ode_rk45`

. Because the states are handled as
`vector`

variables, the solver output is an array of vectors
(`array[] vector`

).

```
transformed parameters {
array[N_t] vector<lower=0>[4] y = ode_rk45(simple_SIR, y0, t0, t,
beta, kappa, gamma, xi, delta);
}
```

The `ode_rk45_tol`

function can be used to specify
tolerances manually:

```
transformed parameters {
array[N_t] vector<lower=0>[4] = ode_rk45_tol(simple_SIR, y0, t0, t,
1e-6, 1e-6, 1000,
beta, kappa, gamma, xi, delta);
}
```

In the old system function, parameters are manually unpacked from
`theta`

and `x_r`

:

```
functions {
array[] real one_comp_mm_elim_abs(real t, array[] real y,
array[] real theta, array[] real x_r,
array[] int x_i) {
array[1] real dydt;
real k_a = theta[1]; // Dosing rate in 1/day
real K_m = theta[2]; // Michaelis-Menten constant in mg/L
real V_m = theta[3]; // Maximum elimination rate in 1/day
real D = x_r[1];
real V = x_r[2];
real dose = 0;
real elim = (V_m / V) * y[1] / (K_m + y[1]);
if (t > 0) {
dose = exp(-k_a * t) * D * k_a / V;
}
dydt[1] = dose - elim;
return dydt;
}
}
```

In the new interface, they are passed directly and so the unpacking is avoided:

```
functions {
vector one_comp_mm_elim_abs(real t,
vector y,
real k_a, // Dosing rate in 1/day
real K_m, // Michaelis-Menten constant in mg/L
real V_m, // Maximum elimination rate in 1/day
real D,
real V) {
vector[1] dydt;
real dose = 0;
real elim = (V_m / V) * y[1] / (K_m + y[1]);
if (t > 0)
dose = exp(- k_a * t) * D * k_a / V;
dydt[1] = dose - elim;
return dydt;
}
}
```

In the old interface the `theta`

and `x_r`

arguments are packed manually, and the `x_i`

argument is
required even though it isn’t used:

```
transformed data {
array[2] real x_r = {D, V};
array[2] int x_i;
}
...
transformed parameters {
array[N_t, 1] real C;
{
array[3] real theta = {k_a, K_m, V_m};
C = integrate_ode_bdf(one_comp_mm_elim_abs, C0, t0, times, theta, x_r, x_i);
}
}
```

In the new interface the arguments are simply passed at the end of
the `ode_bdf`

call (and the `transformed data`

block removed):

```
transformed parameters {
array[N_t] vector[1] mu_C = ode_bdf(one_comp_mm_elim_abs, C0, t0, times,
k_a, K_m, V_m, D, V);
}
```

The `ode_bdf_tol`

function can be used to specify
tolerances manually:

```
transformed parameters {
array[N_t] vector[1] mu_C = ode_bdf_tol(one_comp_mm_elim_abs, C0, t0, times,
1e-8, 1e-8, 1000,
k_a, K_m, V_m, D, V);
}
```