This is an old version, view current version.

13.6 Control Parameters for ODE Solving

For additional control of the solves, both the stiff and non-stiff solvers have function signatures that makes it possible to specify the relative_tolerance, absolute_tolerance, and max_num_steps parameters. These are the same as the regular function names but with _tol appended to the end. All three control arguments must be supplied with this signature (there are no defaults).

vector[2] y_sim[T] = ode_bdf_tol(sho, y0, t0, ts,
                                 relative_tolerance,
                                 absolute_tolerance,
                                 max_num_steps,
                                 theta);

relative_tolerance and absolute_tolerance control the stepsize by specifying a target local error that the solver tries to match, and max_num_steps specifies the maximum number of steps the solver will take between output time points before throwing an error.

The control parameters must be data variables – they cannot be parameters or expressions that depend on parameters, including local variables in any block other than transformed data and generated quantities. User-defined function arguments may be qualified as only allowing data arguments using the data qualifier.

For the RK45 solver, the default values for relative and absolute tolerance are both \(10^{-6}\) and the maximum number of steps between outputs is one million. For the BDF and Adams solvers, the relative and absolute tolerances are \(10^{-10}\) and the maximum number of steps between outputs is one hundred million.

Discontinuous ODE System Function

If there are discontinuities in the ODE system function, it is best to integrate the ODE between the discontinuities, stopping the solver at each one, and restarting it on the other side.

Nonetheless, the ODE solvers will attempt to integrate over discontinuities they encounters in the state function. The accuracy of the solution 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 zero for times \(0 < t < t'\) and then positive for \(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...;

In general it is better to integrate up to t_lag in one solve and then integrate from t_lag onwards in another. Mathematically, the discontinuity can make the problem ill-defined and the numerical integrator may behave erratically around it.

If the location of the discontinuity cannot be controlled precisely, or there is some other rapidly change in ODE behavior, it can be useful to tell the ODE solver to produce output in the neighborhood. This can help the ODE solver avoid indiscriminately stepping over an important feature of the solution.

Tolerance

The relative and absolute tolerance control the accuracy of the solutions generated by the ODE solvers. Relative tolerances are relative to the solution value, whereas absolute tolerances is the maximum absolute error allowed in a solution. Absolute tolerance is more meaningful for ODE solutions that approach zero.

Smaller tolerances produce more accurate solutions, though they require more computation time. The tolerances should be small enough so that setting them lower does not change the statistical properties of posterior samples generated by the Stan program but large enough to avoid unnecessary computation.

Maximum Number of Steps

The maximum number of steps can be used to stop a runaway simulation. This can arise in when MCMC moves to a part of parameter space very far from where a differential equation would typically be solved. In particular this can happen during warmup. With the non-stiff solver, this may happen when the sampler moves to stiff regions of parameter space, which will requires small step sizes.