dopri {dde}  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 = 1e06, atol = 1e06, 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 
func 
Function to integrate. Can be an R function of
arguments 
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 
The output routine; either an R function taking
arguments 
rtol 
The perstep relative tolerance. The total accuracy will be less than this. 
atol 
The perstep 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_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 shortlived fluctuations in
your rhs that the solver may skip over by accident, then specify
a smaller maximum step size here (or use 
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 
step_size_min_allow 
Logical, indicating if when a step size
is driven down to 
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

event_time 
Vector of times to fire events listed in

event_function 
Function to fire at events. For R models
( 
method 
The integration method to use, as a string. The
supported methods are 
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 
verbose 
Be verbose, and print information about each step.
This may be useful for learning about models that misbehave.
Valid values are 
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

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 
return_history 
Logical indicating if history should be
returned alongside the output or discarded. By default, history
is retained if 
dllname 
Name of the shared library (without extension) to
find the function 
parms_are_real 
Logical, indicating if 
ynames 
Logical, indicating if the output should be named
following the names of the input vector 
outnames 
An optional character vector, used when

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 
return_initial 
Logical, indicating if the output should
include the initial conditions. Specifying 
return_time 
Logical, indicating if a row (or column if

return_output_with_y 
Logical, indicating if the output
should be bound together with the returned matrix 
return_statistics 
Logical, indicating if statistics about
the run should be included. If 
restartable 
Logical, indicating if the problem should be
restartable. If 
return_minimal 
Shorthand option  if set to 
obj 
An object to continue from; this must be the results of
running an integration with the option 
copy 
Logical, indicating if the pointer should be copied
before continuing. If 
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 
i 
index within the state vector 
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 510% 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
# 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)