# dopri

##### Integrate ODE/DDE with dopri

Integrate an ODE or DDE with dopri.

##### Usage

```
dopri(y, times, func, parms, ..., n_out = 0L, output = NULL,
rtol = 1e-06, atol = 1e-06, step_size_min = 0,
step_size_max = Inf, step_size_initial = 0, step_max_n = 100000L,
tcrit = NULL, event_time = NULL, event_function = NULL,
method = "dopri5", stiff_check = 0, verbose = FALSE,
callback = NULL, n_history = 0, grow_history = FALSE,
return_history = n_history > 0, dllname = "",
parms_are_real = TRUE, ynames = names(y), outnames = NULL,
return_by_column = TRUE, return_initial = TRUE, return_time = TRUE,
return_output_with_y = TRUE, return_statistics = FALSE,
restartable = FALSE, return_minimal = FALSE)
```dopri5(y, times, func, parms, ...)

dopri853(y, times, func, parms, ...)

dopri_continue(obj, times, y = NULL, ..., copy = FALSE, parms = NULL,
tcrit = NULL, return_history = NULL, return_by_column = NULL,
return_initial = NULL, return_statistics = NULL,
return_time = NULL, return_output_with_y = NULL,
restartable = NULL)

ylag(t, i = NULL)

##### Arguments

- y
Initial conditions for the integration

- times
Times where output is needed. Unlike

`deSolve`

we won't actually stop at these times, but instead interpolate back to get the result.- func
Function to integrate. Can be an R function of arguments

`t, y, parms`

, returning a numeric vector, or it can be the name or address of a C function with arguments`size_t n, double t, const double *y, double *dydt, void *data`

.- parms
Parameters to pass through to the derivatives.

- ...
Dummy arguments - nothing is allowed here, but this means that all further arguments

*must*be specified by name (not order) so I can easily reorder them later on.- n_out
Number of "output" variables (not differential equation variables) to compute via the routine

`output`

.- output
The output routine; either an R function taking arguments

`t, y, parms`

or the name/address of a C function taking arguments`size_t n, double t, const double *y, size_t n_out, double *out, void *data`

.- rtol
The per-step relative tolerance. The total accuracy will be less than this.

- atol
The per-step absolute tolerance.

- step_size_min
The minimum step size. The actual minimum used will be the largest of the absolute value of this

`step_size_min`

or`.Machine$double.eps`

. If the integration attempts to make a step smaller than this, it will throw an error, stopping the integration (note that this differs from the treatment of`hmin`

in`deSolve::lsoda`

).- step_size_max
The largest step size. By default there is no maximum step size (Inf) so the solver can take as large a step as it wants to. If you have short-lived fluctuations in your rhs that the solver may skip over by accident, then specify a smaller maximum step size here (or use

`tcrit`

below).- step_size_initial
The initial step size. By default the integrator will guess the step size automatically, but one can be given here instead.

- step_max_n
The maximum number of steps allowed. If the solver takes more steps than this it will throw an error. Note the number of evaluations of

`func`

will be about 6 times the number of steps (or 11 times if using`method = "dopri853"`

).- tcrit
An optional vector of critical times that the solver must stop at (rather than interpolating over). This can include an end time that we can't go past, or points within the integration that must be stopped at exactly (for example cases where the derivatives change abruptly). Note that this differs from the interpretation of this parameter in deSolve; there

`tcrit`

is a single time that integration may not go past -- with dde we never go past the final time, and this is just for times that fall*within*the range of times in`times`

.- event_time
Vector of times to fire events listed in

`event_function`

at- event_function
Function to fire at events. For R models (

`func`

is an R function and`dllname`

is empty), this must be either a single R function (same function for all events) or a`list`

of R functions. For C models, this must be a singe C function (same requirements as`func`

or`output`

or a list/vector of these as appropriate).- method
The integration method to use, as a string. The supported methods are

`"dopri5"`

(5th order method with 4th order dense output) and`"dopri853"`

(8th order method with 7th order output and embedded 5th and 3rd order schemes). Alternatively, use the functions`dopri5`

or`dopri853`

which simply sets this argument.- stiff_check
How often to check that the problem has become stiff. If zero, then the problem is never checked, and if positive then the problem is checked every

`stiff_check`

accepted steps. The actual check is based off the algorithm in Hairer's implementation of the solvers and may be overly strict, especially for delay equations with the 853 method (in my limited experience with it).- verbose
Be verbose, and print information about each step. This may be useful for learning about models that misbehave. Valid values are

`TRUE`

(enable debugging) or`FALSE`

(disable debugging) or use one of`dopri:::VERBOSE_QUIET`

,`dopri:::VERBOSE_STEP`

or`VERBOSE:::VERBOSE_EVAL`

. If an R function is provided as the argument`callback`

then this function will also be called at each step or evaluation (see below for details).- callback
Callback function that can be used to make verbose output more useful. This can be used to return more information about the evaluation as it proceeds, generally as information printed to the screen. The function must accept arguments

`t`

,`y`

and`dydt`

. See Details for further information.- n_history
Number of history points to retain. This needs to be greater than zero for delay differential equations to work. Alternatively, this may be greater than zero to return model outputs that can be inspected later.

- grow_history
Logical indicating if history should be grown during the simulation. If

`FALSE`

(the default) then when history is used it is overwritten as needed (so only the most recent`n_history`

elements are saved. This may require some tuning so that you have enough history to run your simulation (i.e. to the longest delay) or an error will be thrown when it underflows. The required history length will vary with your delay sizes and with the timestep for dopri. If`TRUE`

, then history will grow as the buffer is exhausted. The growth is geometric, so every time it reaches the end of the buffer it will increase by a factor of about 1.6 (see the`ring`

documentation). This may consume more memory than necessary, but may be useful where you don't want to care about picking the history length carefully.- return_history
Logical indicating if history should be returned alongside the output or discarded. By default, history is retained if

`n_history`

is greater than 0, but that might change (and may not be desirable unless you plan on actually using it).- dllname
Name of the shared library (without extension) to find the function

`func`

(and`output`

if given) in the case where`func`

refers to compiled function.- parms_are_real
Logical, indicating if

`parms`

should be treated as vector of doubles by`func`

(when it is a compiled function). If`TRUE`

(the default), then`REAL(parms)`

, which is`double*`

is passed through. If`FALSE`

then if`params`

is an externalptr type (`EXTPTRSXP`

) we pass through the result of`R_ExternalPtrAddr`

, otherwise we pass`params`

through unmodified as a`SEXP`

. In the last case, in your target function you will need to include`<Rinternals.h>`

, cast to`SEXP`

and then pull it apart using the R API (or Rcpp).- ynames
Logical, indicating if the output should be named following the names of the input vector

`y`

. Alternatively, if`ynames`

is a character vector of the same length as`y`

, these will be used as the output names.- outnames
An optional character vector, used when

`n_out`

is greater than 0, to name the model output matrix.- return_by_column
Logical, indicating if the output should be returned organised by column (rather than row). This incurs a slight cost for transposing the matrices. If you can work with matrices that are transposed relative to

`deSolve`

, then set this to`FALSE`

.- return_initial
Logical, indicating if the output should include the initial conditions. Specifying

`FALSE`

avoids binding this onto the output.- return_time
Logical, indicating if a row (or column if

`return_by_column`

is`TRUE`

) representing time is included. If`FALSE`

, this is not added.- return_output_with_y
Logical, indicating if the output should be bound together with the returned matrix

`y`

(as it is with`deSolve`

). If`FALSE`

, then output will be returned as the attribute`output`

.- return_statistics
Logical, indicating if statistics about the run should be included. If

`TRUE`

, then an integer vector containing the number of target evaluations, steps, accepted steps and rejected steps is returned (the vector is named).- restartable
Logical, indicating if the problem should be restartable. If

`TRUE`

, then the return value of an integration can be passed to`dopri_restart`

to continue the integration after arbitrary changes to the state or the parameters. Note that when using delay models, the integrator is fairly naive about how abrupt changes in the state space are dealt with, and may perform very badly with`method = "dopri853"`

which assumes a fairly smooth problem. Note that this is really only useful for delay differential equations where you want to keep the history but make changes to the parameters or to the state vector while keeping the history of the problem so far.- return_minimal
Shorthand option - if set to

`TRUE`

then it sets all of`return_by_column`

,`return_initial`

,`return_time`

,`return_output_with_y`

to`FALSE`

- obj
An object to continue from; this must be the results of running an integration with the option

`restartable = TRUE`

. Note that continuing a problem moves the pointer along in time (unless`copy = TRUE`

, and that the incoming time (`times[[1]]`

) must equal the previous time*exactly*.- copy
Logical, indicating if the pointer should be copied before continuing. If

`TRUE`

, this is non-destructive with respect to the data in the original pointer so the problem can be restarted multiple times. By default this is`FALSE`

because there is a (potentially very small) cost to this operation.- t
The time to access (not that this is not an offset, but the actual time; within your target function you'd write things like

`tlag(t - 1)`

to get 1 time unit ago.- i
index within the state vector

`y`

to return. The index here is R-style base-1 indexing, so pass`1`

in to access the first element. This can be left`NULL`

to return all the elements or a vector longer than one.

##### Details

Like `deSolve::lsoda`

, this function has *many*
arguments. This is far from ideal, and I would welcome any
approach for simplifying it a bit.

The options `return_by_column`

, `return_initial`

,
`return_time`

, `return_output_with_y`

exist because
these options all carry out modifications of the data at the end
of solving the ODE and this can incur a small but measurable cost.
When solving an ODE repeatedly (e.g., in the context of an MCMC or
optimisation) it may be useful to do as little as possible. For
simple problems this can save around 5-10% of the total
computational time (especially the transpose). The shorthand
option `return_minimal`

will set all to `FALSE`

when
used.

##### Value

At present the return value is transposed relative to deSolve. This might change in future.

##### Verbose output and callbacks

Debugging a failed integration can be difficult, but `dopri`

provides a couple of tools to get more information about where a
failure might have occurred. Most simply, one can pass
`verbose = TRUE`

which will print information about the
time and the step size at each point just before the step is
stated. Passing in `verbose = dde:::VERBOSE_EVAL`

will
print information just before every evaluation of the target
function (there are several evaluations per step).

However, this does not provide information about the state just
before failure. To get that, one must provide a `callback`

function - this is an R function that will be called just before
a step or evaluation (based on the value of the `verbose`

argument) in place of the default print. Define a callback
function with arguments `t`

, `h`

and `y`

where
`t`

is the time (beginning of a step or location of an
evaluation), `h`

is the step size (or `NA`

for an
evaluation) and `y`

is the state at the point of the step
or evaluation. Your callback function can do anything - you can
print to the screen (using `cat`

or `message`

), you
can store results using a closure and `<<-`

or you could
conditionally use a `browser()`

call to debug
interactively. However, it is not possible for the callback to
affect the solution short of throwing an error and interrupting
it. See the Examples for an example of use.

##### See Also

`dopri_interpolate`

which can be used to
efficiently sample from output of `dopri`

, and the package
vignette which shows in more detail how to solve delay
differential equations and to use compiled objective functions.

##### Examples

```
# NOT RUN {
# The lorenz attractor:
lorenz <- function(t, y, p) {
sigma <- p[[1L]]
R <- p[[2L]]
b <- p[[3L]]
c(sigma * (y[[2L]] - y[[1L]]),
R * y[[1L]] - y[[2L]] - y[[1L]] * y[[3L]],
-b * y[[3L]] + y[[1L]] * y[[2L]])
}
p <- c(10, 28, 8 / 3)
y0 <- c(10, 1, 1)
tt <- seq(0, 100, length.out = 40000)
y <- dde::dopri(y0, tt, lorenz, p, return_time = FALSE)
plot(y[, c(1, 3)], type = "l", lwd = 0.5, col = "#00000066")
# If we want to print progress as the integration progresses we can
# use the verbose argument:
y <- dde::dopri(y0, c(0, 0.1), lorenz, p, verbose = TRUE)
# Or print the y values too using a callback:
callback <- function(t, h, y) {
message(sprintf("t: %f, h: %e, y: [%s]", t, h,
paste(format(y, 5), collapse = ", ")))
}
y <- dde::dopri(y0, c(0, 0.1), lorenz, p, verbose = TRUE,
callback = callback)
# }
```

*Documentation reproduced from package dde, version 1.0.0, License: MIT + file LICENSE*