13.6 Control parameters for ODE solving
For additional control of the solves, both the stiff and non-stiff
forward ODE 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).
array[T] vector[2] y_sim = ode_bdf_tol(sho, y0, t0, ts,
relative_tolerance,
absolute_tolerance,
max_num_steps,
theta);relative_tolerance and absolute_tolerance control accuracy the solver tries to achieve, 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 and Cash-Karp solvers, 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, 0]';
} else {
// ... return non-zero vector...
}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 tolerance RTOL and absolute tolerance ATOL control the accuracy of the numerical solution. Specifically, when solving an ODE with unknowns \(y=(y_1,\dots,y_n)^T\), at every step the solver controls estimated local error \(e=(e_1,\dots,e_n)^T\) through its weighted root-mean-square norm (Serban and Hindmarsh (2005), Hairer, Nørsett, and Wanner (1993))
\[ \sqrt{\sum_{i=1}^n{\frac{1}{n}\frac{e_i^2}{(\text{RTOL}\times y_i + \text{ATOL})^2}}} < 1 \] by reducing the stepsize when the inequality is not satisfied.
To understand the roles of the two tolerances it helps to assume \(y\) at
opposite scales in the above expression: on one hand the absolute
tolerance has little effect when \(y_i \gg 1\), on the other the
relative tolerance can
not affect the norm when \(y_i = 0\). Users are strongly encouraged to carefully choose
tolerance values according to the ODE and its application. One can follow
Brenan, Campbell, and Petzold (1995) for a rule of thumb:
let \(m\) be the number of significant digits required for \(y\), set
\(\text{RTOL}=10^{-(m+1)}\), and set ATOL at
which \(y\) becomes insignificant. Note that the same weighted root-mean-square norm
is used to control nonlinear solver convergence in bdf and adams solvers, and the same
tolerances are used to control forward sensitivity calculation. See
Serban and Hindmarsh (2005) for details.
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.