## 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 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 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 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 toleranes 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.

*References*

Brenan, K. E., S. L. Campbell, and L. R. Petzold. 1995. *Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations*. SIAM Classics in Applied Mathematics 14. Philadelphia: Society for Industrial; Applied Mathematics.

Hairer, Ernst, Syvert P. Nørsett, and Gerhard Wanner. 1993. *Solving Ordinary Differential Equations I: Nonstiff Problems*. 2nd ed. Springer Series in Computational Mathematics, Springer Ser.Comp.Mathem. Hairer,E.:Solving Ordinary Diff. Berlin Heidelberg: Springer-Verlag.

Serban, Radu, and Alan C Hindmarsh. 2005. “CVODES: The Sensitivity-Enabled ODE Solver in SUNDIALS.” In *ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference*, 257–69. American Society of Mechanical Engineers.