dopri: Integrate ODE/DDE with dopri
In dde: Solve Delay Differential Equations

dopri

R Documentation

Integrate ODE/DDE with dopri

Description

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,
  step_size_min_allow = FALSE, 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 by default, stopping the integration (note that this differs from the treatment of `hmin` in `deSolve::lsoda`). See `allow_step_size_min` to change this behaviour.
`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"`).
`step_size_min_allow`	Logical, indicating if when a step size is driven down to `step_size_min` we should allow it to proceed. This is the behaviour in of `hmin` in `deSolve::lsoda`.
`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.

Examples


# 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)

dde documentation built on Sept. 23, 2024, 5:09 p.m.