| Title: | ODE Solver with Automatic Differentiation, in C++ Header Files |
|---|---|
| Description: | odelia provides an ODE solver in C++ header files, offering both an explicit adaptive-step Runge-Kutta 4-5 method and an implicit, adaptive-step RODAS4(3) Rosenbrock method for stiff systems, with an interface to R via Rcpp. ODE systems can be templated on their scalar type to support automatic differentiation, enabling exact reverse-mode gradients of a recorded run and exact Jacobians for the implicit solver. A system written as an R function can be solved by the same steppers. |
| Authors: | Daniel Falster [aut, cre] (ORCID: <https://orcid.org/0000-0002-9814-092X>), Richard FitzJohn [aut], Andrew O'Reilly-Nugent [aut] |
| Maintainer: | Daniel Falster <[email protected]> |
| License: | AGPL (>= 3) |
| Version: | 0.11.1 |
| Built: | 2026-11-13 01:19:39 UTC |
| Source: | https://github.com/traitecoevo/odelia |
odelia provides an ODE solver in C++ header files, offering both an explicit adaptive-step Runge-Kutta 4-5 method and an implicit, adaptive-step RODAS4(3) Rosenbrock method for stiff systems, with an interface to R via Rcpp. ODE systems can be templated on their scalar type to support automatic differentiation, enabling exact reverse-mode gradients of a recorded run and exact Jacobians for the implicit solver. A system written as an R function can be solved by the same steppers.
The DESCRIPTION file:
| Package: | odelia |
| Type: | Package |
| Title: | ODE Solver with Automatic Differentiation, in C++ Header Files |
| Version: | 0.11.1 |
| Authors@R: | c( person("Daniel", "Falster", , "[email protected]", c("aut", "cre"), comment = c(ORCID = "0000-0002-9814-092X")), person("Richard", "FitzJohn", , "", "aut"), person("Andrew", "O'Reilly-Nugent", , "", "aut") ) |
| Maintainer: | Daniel Falster <[email protected]> |
| Description: | odelia provides an ODE solver in C++ header files, offering both an explicit adaptive-step Runge-Kutta 4-5 method and an implicit, adaptive-step RODAS4(3) Rosenbrock method for stiff systems, with an interface to R via Rcpp. ODE systems can be templated on their scalar type to support automatic differentiation, enabling exact reverse-mode gradients of a recorded run and exact Jacobians for the implicit solver. A system written as an R function can be solved by the same steppers. |
| Depends: | R (>= 4.2.0) |
| Imports: | Rcpp (>= 1.0.10), R6 (>= 2.5.1), dplyr (>= 1.1.0), tibble (>= 3.1.0) |
| Suggests: | deSolve, here, knitr, rmarkdown, rootSolve, testthat (>= 3.0.0), withr |
| LinkingTo: | Rcpp |
| VignetteBuilder: | knitr |
| Encoding: | UTF-8 |
| License: | AGPL (>= 3) |
| URL: | https://github.com/traitecoevo/odelia, https://traitecoevo.github.io/odelia/ |
| BugReports: | https://github.com/traitecoevo/odelia/issues |
| Config/testthat/edition: | 3 |
| Config/roxygen2/version: | 8.0.0 |
| Repository: | https://traitecoevo.r-universe.dev |
| Date/Publication: | 2026-11-12 22:23:50 UTC |
| RemoteUrl: | https://github.com/traitecoevo/odelia |
| RemoteRef: | HEAD |
| RemoteSha: | ed4fc29dbe05c73957767a563efda5730239d248 |
| Author: | Daniel Falster [aut, cre] (ORCID: <https://orcid.org/0000-0002-9814-092X>), Richard FitzJohn [aut], Andrew O'Reilly-Nugent [aut] |
Index of help topics:
domain_error Signal that a state is outside the model's
domain
Drivers Drivers R6 Class
Lorenz_Solver Lorenz Solver R6 Class
LorenzSystem Lorenz System R6 Class
ode_counts What a solve cost
ode_solve Solve an ODE given as an R function, deSolve
style
ode_steady_state Equilibrium of an ODE given as an R function,
and how it moves with the parameters
OdeControl OdeControl R6 Class
odelia_load_dll Load odelia shared library with configurable
symbol visibility
odelia-package ODE Solver with Automatic Differentiation, in
C++ Header Files
OdeSolver An ODE solver over an R right-hand side
Daniel Falster <[email protected]>
Daniel Falster [aut, cre] (ORCID: <https://orcid.org/0000-0002-9814-092X>), Richard FitzJohn [aut], Andrew O'Reilly-Nugent [aut]
Call this from a right-hand side, Jacobian or validity predicate handed to [OdeSolver] or [ode_solve()], at any depth below it, to say that the state it was given is one the model has no meaning for – a pool below zero, a probability above one, an inner solve that did not converge – as distinct from an error in the code. The callback is left at once, as by 'return()' (its 'on.exit()' expressions run), and the solver *rejects the current step* and retries it smaller, exactly as it does for a compiled system throwing 'util::DomainError'. Any other error raised in a callback propagates unchanged and ends the solve, which is what a bug should do.
domain_error(message)domain_error(message)
message |
What is wrong with the state. |
If the smallest permitted step still leaves the domain, the solve stops with an error naming this message, so make it say what was wrong. Called outside a solver callback, this is an ordinary error of class 'odelia_domain_error'.
Does not return.
rhs <- function(t, y) { if (y[1] < 0) domain_error("y must stay non-negative") -sqrt(y[1]) }rhs <- function(t, y) { if (y[1] < 0) domain_error("y must stay non-negative") -sqrt(y[1]) }
R6 wrapper for external forcing data (drivers)
ptrExternal pointer to the underlying C++ Drivers object.
Drivers$new()Initialize a 'Drivers' instance.
Drivers$new()
Drivers$set_constant()Set a constant driver value.
Drivers$set_constant(name, value)
nameDriver variable name.
valueConstant driver value.
Drivers$set_variable()Set a variable driver from paired vectors.
Drivers$set_variable(name, x, y)
nameDriver variable name.
xNumeric vector of input coordinates (e.g., time).
yNumeric vector of values corresponding to 'x'.
Drivers$set_extrapolate()Configure extrapolation behavior for a driver.
Drivers$set_extrapolate(name, extrapolate)
nameDriver variable name.
extrapolateLogical flag controlling extrapolation beyond data range.
Drivers$evaluate()Evaluate a driver at one input value.
Drivers$evaluate(name, x)
nameDriver variable name.
xNumeric vector of input coordinates (e.g., time).
Drivers$evaluate_range()Evaluate a driver over a vector of input values.
Drivers$evaluate_range(name, x)
nameDriver variable name.
xNumeric vector of input coordinates (e.g., time).
Drivers$get_names()Return names of configured drivers.
Drivers$get_names()
Drivers$clear()Remove all configured drivers.
Drivers$clear()
Drivers$clone()The objects of this class are cloneable with this method.
Drivers$clone(deep = FALSE)
deepWhether to make a deep clone.
R6 wrapper for the Lorenz ODE solver.
ptrExternal pointer to the underlying C++ solver object.
Lorenz_Solver$new()Initialize a solver for a Lorenz system. The solver takes
its own copy of the system as it is at this moment and steps that copy:
set_params() or set_state() on the LorenzSystem
afterwards do not reach the solver. Use the solver's own
$set_params() and $set_state().
Lorenz_Solver$new(System_xp, control_xp, method = "rkck")
System_xpExternal pointer to the Lorenz system object.
control_xpExternal pointer to the ODE control object.
methodIntegration method: "rkck" (default, explicit Cash-Karp
RK 4(5)), "dopri" (explicit Dormand-Prince 5(4)) or "rodas"
(implicit RODAS4(3) Rosenbrock, for stiff systems).
Lorenz_Solver$pars()Parameters of the system this solver steps.
Lorenz_Solver$pars()
Lorenz_Solver$set_params()Set the parameters this solver steps with, and re-seed it
at its current state and time as $set_state() would. Call
$reset() as well to start again from the initial state.
Lorenz_Solver$set_params(params)
paramsNumeric vector of system parameters, one entry per parameter in the system's order (for Lorenz: sigma, R, b).
Lorenz_Solver$time()Get current solver time.
Lorenz_Solver$time()
Lorenz_Solver$state()Get current solver state.
Lorenz_Solver$state()
Lorenz_Solver$set_state()Set current solver state and time.
Lorenz_Solver$set_state(y, time)
yNumeric state vector.
timeScalar time value.
Lorenz_Solver$times()Get stored solver times.
Lorenz_Solver$times()
Lorenz_Solver$advance_adaptive()Advance solver using adaptive stepping.
Lorenz_Solver$advance_adaptive(times)
timesNumeric vector of requested output times.
Lorenz_Solver$advance_fixed()Advance solver using fixed stepping.
Lorenz_Solver$advance_fixed(times)
timesNumeric vector of requested output times.
Lorenz_Solver$advance_euler()Advance solver using fixed-step forward Euler.
Lorenz_Solver$advance_euler(times)
timesNumeric vector of requested output times.
Lorenz_Solver$step()Advance solver by one step.
Lorenz_Solver$step()
Lorenz_Solver$reset()Reset solver to its initial state.
Lorenz_Solver$reset()
Lorenz_Solver$collect()Get or set history collection behavior.
Lorenz_Solver$collect(value)
valueLogical flag for collect-history behavior.
Lorenz_Solver$history_size()Return number of stored history entries.
Lorenz_Solver$history_size()
Lorenz_Solver$history_step()Return one history entry by index.
Lorenz_Solver$history_step(i)
iInteger index into solver history.
Lorenz_Solver$history()Return history as a tibble.
Lorenz_Solver$history()
Lorenz_Solver$set_target()Set the calibration target for
$loss_and_gradient().
Lorenz_Solver$set_target(schedule, target, obs_indices)
scheduleNumeric vector: the step times
$loss_and_gradient() replays, one entry per step. Normally a
reference run's $times(), so each evaluation takes the same steps
and the loss is a smooth function of the inputs.
targetNumeric matrix of observed states, one row per entry of
obs_indices and one column per state variable.
obs_indicesInteger indices (1-based) into schedule at which
each row of target was observed.
Lorenz_Solver$loss_and_gradient()Least-squares loss against the target, and its exact
gradient by reverse-mode automatic differentiation. The schedule is
replayed under this solver's method, except that a
"rodas" solver's is replayed under "rkck": the
Rosenbrock stepper cannot be swept.
Lorenz_Solver$loss_and_gradient(initial_state = NULL, params = NULL)
initial_stateInitial state to evaluate at, or NULL to keep
the system's.
paramsParameters to evaluate at, or NULL to keep the
system's.
A list with loss (sum of squared differences) and
gradient: d(loss)/d(params), then d(loss)/d(initial_state), for
whichever were given.
Lorenz_Solver$clone()The objects of this class are cloneable with this method.
Lorenz_Solver$clone(deep = FALSE)
deepWhether to make a deep clone.
R6 wrapper for Lorenz system
ptrExternal pointer to the underlying C++ Lorenz system object.
LorenzSystem$new()Initialize a Lorenz system object.
LorenzSystem$new(sigma, R, b)
sigmaLorenz parameter sigma.
RLorenz parameter R.
bLorenz parameter b.
LorenzSystem$pars()Return current system parameters.
LorenzSystem$pars()
LorenzSystem$set_params()Set model parameters.
LorenzSystem$set_params(params)
paramsNumeric vector of system parameters.
LorenzSystem$set_state()Set system state and time.
LorenzSystem$set_state(y, time = 0)
yNumeric state vector.
timeScalar time value.
LorenzSystem$state()Return current state.
LorenzSystem$state()
LorenzSystem$rates()Return current rates.
LorenzSystem$rates()
LorenzSystem$set_initial_state()Set initial state and initial time.
LorenzSystem$set_initial_state(y, t0 = 0)
yNumeric state vector.
t0Initial time value.
LorenzSystem$reset()Reset the system to its initial condition.
LorenzSystem$reset()
LorenzSystem$clone()The objects of this class are cloneable with this method.
LorenzSystem$clone(deep = FALSE)
deepWhether to make a deep clone.
The number of right-hand-side evaluations ('n_rhs'), Jacobian formations ('n_jac'), accepted steps ('n_steps') and rejected attempts ('n_rejections') behind a result of [ode_solve()], as a named vector. For a right-hand side that is itself expensive these are the whole cost of the integration; 'n_rejections' is what the step-size controller wasted.
ode_counts(x, ...) ## S3 method for class 'odelia_solution' ode_counts(x, ...) ## S3 method for class 'odelia_steady_state' ode_counts(x, ...)ode_counts(x, ...) ## S3 method for class 'odelia_solution' ode_counts(x, ...) ## S3 method for class 'odelia_steady_state' ode_counts(x, ...)
x |
A result of [ode_solve()]. |
... |
Ignored. |
A named numeric vector.
out <- ode_solve(function(t, y, p) -y, y0 = 1, times = c(0, 1)) ode_counts(out)out <- ode_solve(function(t, y, p) -y, y0 = 1, times = c(0, 1)) ode_counts(out)
A convenience over [OdeSolver] shaped like 'deSolve::ode()': the right-hand side is 'func(t, y, parms)' returning a list whose first element is the rates (a bare numeric vector is accepted too), the optional Jacobian is 'jacfunc(t, y, parms)' returning an n by n matrix with column j = 'd(dy/dt)/dy[j]', and the result is a matrix with a 'time' column and one column per state variable, one row per requested time. A function written for deSolve runs unchanged.
ode_solve( func, y0, times, parms = NULL, jacfunc = NULL, method = "dopri", control = NULL, rtol = 1e-06, atol = 1e-06, controller = "hairer", autonomous = FALSE, jac_fd_step = 1e-06, jac_fd_floor = 1e-05, dense = NULL )ode_solve( func, y0, times, parms = NULL, jacfunc = NULL, method = "dopri", control = NULL, rtol = 1e-06, atol = 1e-06, controller = "hairer", autonomous = FALSE, jac_fd_step = 1e-06, jac_fd_floor = 1e-05, dense = NULL )
func |
The right-hand side, 'func(t, y, parms)'. |
y0 |
Initial state; names, if any, become column names. |
times |
Times to report at; the first is the initial time. |
parms |
Passed through to 'func' and 'jacfunc'. |
jacfunc |
'NULL', or the Jacobian 'jacfunc(t, y, parms)'. |
method |
'"dopri"' (explicit Dormand–Prince 5(4), the default, as ‘deSolve'’s 'ode45'), '"rkck"' (explicit Cash–Karp 4(5)) or '"rodas"' (implicit RODAS4(3), for stiff problems). |
control |
An [OdeControl], or 'NULL' for the defaults with the tolerances below, the 'controller' below, and the largest step set to the span of 'times' (an [OdeControl] made directly caps the step at 10). With a 'control' given, 'rtol', 'atol' and 'controller' are ignored. |
rtol, atol
|
Shorthand for the control's relative and absolute tolerances when 'control' is 'NULL'. |
controller |
The step-size rule when 'control' is 'NULL': '"hairer"' (the default) or '"gsl"'; see [OdeControl]. |
autonomous |
'TRUE' if 'func' does not depend on 't'. |
jac_fd_step, jac_fd_floor
|
The finite-difference Jacobian's relative step and small-component floor; see [OdeSolver]. |
dense |
'TRUE' to let the stepper choose its own steps and read the requested times off the interpolant of the step spanning each, as ‘deSolve'’s 'lsoda' does; 'FALSE' to land a step on every requested time, as its 'ode45' does, which costs more steps when the output is finer than the steps. 'NULL', the default, means 'TRUE' under '"dopri"', whose interpolant has the stepper's order, and 'FALSE' under '"rkck"' and '"rodas"', whose interpolant is cubic Hermite, one order short (on Lorenz at 1e-6 the Cash–Karp dense rows err by 6e-4 where landed rows err by 3e-11, measured). |
A numeric matrix of class 'odelia_solution', 'length(times)' rows; [ode_counts()] on it says what the solve cost.
lorenz <- function(t, y, p) { list(c(p[["sigma"]] * (y[2] - y[1]), y[1] * (p[["rho"]] - y[3]) - y[2], y[1] * y[2] - p[["beta"]] * y[3])) } out <- ode_solve(lorenz, y0 = c(x = 1, y = 1, z = 1), times = seq(0, 2, by = 0.1), parms = c(sigma = 10, rho = 28, beta = 8 / 3), autonomous = TRUE) head(out) ode_counts(out)lorenz <- function(t, y, p) { list(c(p[["sigma"]] * (y[2] - y[1]), y[1] * (p[["rho"]] - y[3]) - y[2], y[1] * y[2] - p[["beta"]] * y[3])) } out <- ode_solve(lorenz, y0 = c(x = 1, y = 1, z = 1), times = seq(0, 2, by = 0.1), parms = c(sigma = 10, rho = 28, beta = 8 / 3), autonomous = TRUE) head(out) ode_counts(out)
Solve 'func(t, y, parms) = 0' for a fixed point 'y*' by damped Newton, say whether it is attracting, and take its sensitivity to every numeric parameter by the implicit function theorem, 'dy*/dp = -(df/dy)^-1 (df/dp)', with no integration through the transient. The right-hand side has the shape [ode_solve()] takes, so a function written for 'deSolve' (or for 'rootSolve::steady()') runs unchanged.
ode_steady_state( func, y0, parms = NULL, jacfunc = NULL, t0 = 0, tol = 1e-10, max_iter = 100L, warmup = NULL, method = "rodas", control = NULL, rtol = 1e-06, atol = 1e-06, sensitivity = TRUE, parms_fd_step = 1e-06, autonomous = FALSE, jac_fd_step = 1e-06, jac_fd_floor = 1e-05, line_search = TRUE )ode_steady_state( func, y0, parms = NULL, jacfunc = NULL, t0 = 0, tol = 1e-10, max_iter = 100L, warmup = NULL, method = "rodas", control = NULL, rtol = 1e-06, atol = 1e-06, sensitivity = TRUE, parms_fd_step = 1e-06, autonomous = FALSE, jac_fd_step = 1e-06, jac_fd_floor = 1e-05, line_search = TRUE )
func |
The right-hand side, 'func(t, y, parms)', returning the rates or a list whose first element is the rates. |
y0 |
Starting guess; names, if any, name the result. |
parms |
Passed through to 'func' and 'jacfunc'. Its numeric elements are what the sensitivity is taken with respect to. |
jacfunc |
'NULL', or the Jacobian 'jacfunc(t, y, parms)'. |
t0 |
The time 'func' is evaluated at. |
tol |
Convergence: the largest absolute rate at the solution. |
max_iter |
The most Newton iterations to take. |
warmup |
'NULL' for Newton alone, or the span of time to integrate the transient over before retrying Newton when its root is not attracting; a vector of at least two times, the first 't0', is taken as the times to land on instead. |
method |
The stepper for the warm-up: '"rodas"' (implicit RODAS4(3), the default, for a stiff system), '"dopri"' or '"rkck"'. |
control |
An [OdeControl] for the warm-up, or 'NULL' for the defaults with 'rtol' and 'atol'. |
rtol, atol
|
The warm-up's tolerances when 'control' is 'NULL'. |
sensitivity |
'TRUE' for every numeric parameter, 'FALSE' for none, or the names or indices of the parameters to take it for. |
parms_fd_step |
Relative step of the central difference in each parameter: 'h = parms_fd_step * max(abs(p), 1)'. |
autonomous |
'TRUE' if 'func' does not depend on 't', which saves the evaluation behind 'time_dependence'. |
jac_fd_step, jac_fd_floor
|
The finite-difference Jacobian's relative step and small-component floor; see [OdeSolver]. |
line_search |
'TRUE' to backtrack on the residual when a full Newton step would not reduce it, which widens the basin of convergence. |
Newton needs the Jacobian 'df/dy': supply 'jacfunc(t, y, parms)' returning the n by n matrix whose column j is the derivative with respect to 'y[j]', or leave it 'NULL' and it is formed by forward differences at 'n' right-hand-side evaluations per iteration. Each iteration also costs one evaluation per line-search trial, usually one.
Newton finds *a* root, not necessarily an attracting one: from a guess near the trivial equilibrium of a demographic model it converges there in a step or two, and 'stable' says so. Give 'warmup' to seek an attractor instead: when Newton's root is not attracting, or Newton does not converge, the transient is integrated from 'y0' for that span with 'method' (RODAS for a stiff system) and Newton is retried from where it ends; 'warmed' reports that this happened.
The sensitivity is to the parameters 'func' can be perturbed in from R: every element of a numeric 'parms', or every element of a list 'parms' that is a single number. Each column is a central difference of 'func' at 'y*', two evaluations per parameter, so its accuracy is that of the difference ('parms_fd_step'), not of the Jacobian. 'df/dy' is singular at a bifurcation point, where the theorem does not apply; the sensitivity stops with an error there.
Equilibrium is only meaningful for an autonomous system. 'time_dependence' is the size of 'df/dt' at the solution, formed by one extra evaluation unless 'autonomous = TRUE'; a value above finite-difference noise means 'func' depends on 't' and the "equilibrium" is a root at 't0' only.
A list of class 'odelia_steady_state': 'y', the equilibrium estimate (named as 'y0'); 'residual', the rates there, and 'residual_norm', the largest absolute one; 'iterations', 'converged' and 'warmed'; 'jacobian', the n by n 'df/dy' at 'y*', column j the derivative with respect to 'y[j]'; 'eigenvalues', those of the Jacobian, complex; 'spectral_abscissa', the largest real part; and 'stable', 'TRUE' when that is negative; 'time_dependence', the largest absolute 'df/dt' at the solution; 'sensitivity', the n by p matrix 'dy*/dp', rows named as 'y', columns as the parameters, or 'NULL' when none was asked for or possible; and 'counts', what the solve cost ([ode_counts()]). Without convergence the Jacobian and everything read off it are 'NULL', and a warning says so.
# y0' = a - b y0 ; y1' = y0^2 - c y1 : the fixed point is (a/b, a^2/(b^2 c)). f <- function(t, y, p) c(p[["a"]] - p[["b"]] * y[1], y[1]^2 - p[["c"]] * y[2]) eq <- ode_steady_state(f, y0 = c(n = 0, m = 0), parms = c(a = 2, b = 1.5, c = 0.7)) eq$y eq$stable eq$sensitivity # dy*/da, dy*/db, dy*/dc # A logistic population has a trivial root at 0 that Newton finds from a # small guess; warmup integrates past it to the attractor. g <- function(t, y, p) p[["r"]] * y * (1 - y / p[["K"]]) ode_steady_state(g, y0 = 1e-3, parms = c(r = 1, K = 10))$y ode_steady_state(g, y0 = 1e-3, parms = c(r = 1, K = 10), warmup = 50)$y# y0' = a - b y0 ; y1' = y0^2 - c y1 : the fixed point is (a/b, a^2/(b^2 c)). f <- function(t, y, p) c(p[["a"]] - p[["b"]] * y[1], y[1]^2 - p[["c"]] * y[2]) eq <- ode_steady_state(f, y0 = c(n = 0, m = 0), parms = c(a = 2, b = 1.5, c = 0.7)) eq$y eq$stable eq$sensitivity # dy*/da, dy*/db, dy*/dc # A logistic population has a trivial root at 0 that Newton finds from a # small guess; warmup integrates past it to the attractor. g <- function(t, y, p) p[["r"]] * y * (1 - y / p[["K"]]) ode_steady_state(g, y0 = 1e-3, parms = c(r = 1, K = 10))$y ode_steady_state(g, y0 = 1e-3, parms = c(r = 1, K = 10), warmup = 50)$y
R6 wrapper for ODE solver control settings
ptrExternal pointer to the underlying C++ OdeControl object.
OdeControl$new()Initialize an 'OdeControl' instance.
OdeControl$new()
OdeControl$get_controls()Get current solver control settings.
OdeControl$get_controls()
OdeControl$set_controls()Set all solver control settings at once.
OdeControl$set_controls( tol_abs, tol_rel, a_y, a_dydt, step_size_min, step_size_max, step_size_initial )
tol_absAbsolute tolerance.
tol_relRelative tolerance.
a_yState scaling coefficient.
a_dydtDerivative scaling coefficient.
step_size_minMinimum allowed step size.
step_size_maxMaximum allowed step size.
step_size_initialInitial step size.
OdeControl$set_tol_abs()Set absolute tolerance.
OdeControl$set_tol_abs(tol_abs)
tol_absAbsolute tolerance.
OdeControl$set_tol_rel()Set relative tolerance.
OdeControl$set_tol_rel(tol_rel)
tol_relRelative tolerance.
OdeControl$set_a_y()Set state scaling coefficient.
OdeControl$set_a_y(a_y)
a_yState scaling coefficient.
OdeControl$set_a_dydt()Set derivative scaling coefficient.
OdeControl$set_a_dydt(a_dydt)
a_dydtDerivative scaling coefficient.
OdeControl$set_step_size_min()Set minimum step size.
OdeControl$set_step_size_min(step_size_min)
step_size_minMinimum allowed step size.
OdeControl$set_step_size_max()Set maximum step size.
OdeControl$set_step_size_max(step_size_max)
step_size_maxMaximum allowed step size.
OdeControl$set_controller()Set the step-size rule, '"gsl"' or '"hairer"'.
OdeControl$set_controller(controller)
controllerWhich rule turns the error estimate into the next step: '"hairer"' (the default; the error in the RMS norm over components, accept at 1 and rescale after every step, as Hairer's codes and ‘deSolve' do, damped by the previous step’s error as 'dopri5.f' is and ‘deSolve'’s 'rk' is not; see 'set_stabilisation()') or '"gsl"' (the worst component's error, accept below 1.1, grow only below 0.5, otherwise keep the step). '"gsl"' rejects many more attempts; ‘"hairer"'’s RMS norm bounds each component more loosely when the state is large and the error sits in a few components.
OdeControl$get_controller()The step-size rule in force.
OdeControl$get_controller()
OdeControl$set_stabilisation()Set the ‘"hairer"' rule’s stabilisation, 'beta' in 'dopri5.f': 0.04 (the default) damps step growth while the error is rising, which keeps an explicit stepper off its stability limit; 0 is the classical rule, which takes fewer steps on smooth problems.
OdeControl$set_stabilisation(stabilisation)
stabilisationA number in '[0, 0.2)'.
OdeControl$get_stabilisation()The ‘"hairer"' rule’s stabilisation in force.
OdeControl$get_stabilisation()
OdeControl$set_step_size_initial()Set initial step size.
OdeControl$set_step_size_initial(step_size_initial)
step_size_initialInitial step size.
OdeControl$clone()The objects of this class are cloneable with this method.
OdeControl$clone(deep = FALSE)
deepWhether to make a deep clone.
Locates the installed 'odelia' shared library and calls [base::dyn.load()] on it. This is useful in 'Rcpp::sourceCpp()' workflows where symbols from the package shared object must be available to a temporary shared object.
odelia_load_dll(local = FALSE, now = TRUE)odelia_load_dll(local = FALSE, now = TRUE)
local |
Passed to [base::dyn.load()]. Default 'FALSE' to expose symbols globally for downstream dynamic linking. |
now |
Passed to [base::dyn.load()]. Default 'TRUE'. |
Invisibly returns the resolved path to the 'odelia' shared library.
Integrate a system whose right-hand side is an R function, with odelia's adaptive step control and any of its steppers: the explicit Dormand–Prince 5(4) pair (the method behind ‘deSolve'’s 'ode45', with dense output of its own order), the explicit Cash–Karp 4(5) pair, or the implicit RODAS4(3) Rosenbrock stepper for stiff problems. The solver is driven a step at a time (or to a sequence of times), and the state can be reset between steps – at a different length if need be – which is what a consumer with events needs.
The right-hand side 'rhs(t, y)' returns the rates as a numeric vector the length of 'y'; with 'parms' given it is called as 'rhs(t, y, parms)', and a list whose first element is the rates is accepted, so a function written for 'deSolve' runs unchanged. The implicit stepper needs the Jacobian 'd(dy/dt)/dy': supply 'jac(t, y)' (or 'jac(t, y, parms)') returning the n by n matrix whose column j is the derivative with respect to 'y[j]', or leave it 'NULL' and it is formed by forward differences at 'n' extra right-hand-side evaluations per step.
Only the number of right-hand-side evaluations matters for a right-hand side that is itself expensive, and '$counts()' reports it. Per accepted step RODAS makes six (five stages and the derivative at the new point), plus one for the time derivative unless 'autonomous = TRUE', plus the Jacobian once – 'n' evaluations by finite differences, one call to 'jac' otherwise – which is kept across a retry of a rejected step. A rejected attempt costs its six stage evaluations. For a cheap right-hand side the cost is instead one R call per stage, which is why an R right-hand side runs no faster here than through 'deSolve'; see the Lorenz benchmark.
A callback that finds its state impossible calls [domain_error()] and the step is rejected and retried smaller. Any other error propagates and leaves the solver holding a half-finished step; 'set_state()' makes it usable again. After every accepted 'step()', and after 'set_state()', the last call to 'rhs' was at exactly '(time(), state())', so a closure that keeps extra results from its last evaluation can rely on them matching the solver's state.
OdeSolver$new()Create a solver over an R right-hand side. Evaluates 'rhs' once, at '(t0, y0)'.
OdeSolver$new( rhs, y0, t0 = 0, jac = NULL, state_valid = NULL, parms = NULL, control = NULL, method = "dopri", autonomous = FALSE, jac_fd_step = 1e-06, jac_fd_floor = 1e-05 )
rhsFunction 'rhs(t, y)' returning 'dy/dt', a numeric vector the length of 'y'.
y0Initial state (numeric, length at least one).
t0Initial time.
jac'NULL', or a function 'jac(t, y)' returning the n by n Jacobian matrix with column j = 'd(dy/dt)/dy[j]'.
state_valid'NULL', or a predicate 'state_valid(t, y)' returning 'TRUE' for a state the model accepts; a step landing on a refused state is rejected.
parms'NULL', or a value passed as a third argument to every callback.
controlAn [OdeControl], or 'NULL' for the defaults.
method'"dopri"' (explicit Dormand–Prince 5(4), the default; the one with dense output of its own order), '"rkck"' (explicit Cash–Karp 4(5)) or '"rodas"' (implicit RODAS4(3), for stiff problems).
autonomous'TRUE' if 'rhs' does not depend on 't', which saves the implicit stepper one evaluation per step for the time derivative.
jac_fd_stepRelative step for the finite-difference Jacobian: 'h_j = jac_fd_step * max(abs(y[j]), jac_fd_floor)'. A right-hand side that is itself an iterative solve wants a larger step than the default to stay above its own noise.
jac_fd_floorThe size below which a state component counts as small for that perturbation (the default 1e-5 is ‘rodas.f'’s). A floor of 1 perturbs a component of size 1e-5 by a tenth of itself; a floor as small as a tight absolute tolerance makes the perturbation vanish in the subtraction. Lower it only for a state whose components legitimately live below 1e-5. A component sitting on the upper edge of its domain is perturbed downwards instead when the right-hand side refuses the upward point with [domain_error()].
OdeSolver$step()Take one adaptive step, not passing 'time_max'. Refused when already at a finite 'time_max'.
OdeSolver$step(time_max = Inf)
time_maxA time the step must not pass; 'Inf' for no bound.
OdeSolver$advance_adaptive()Advance to each of 'times' in turn by adaptive steps, landing on each exactly.
OdeSolver$advance_adaptive(times)
timesTimes to advance to; the first must be the current time.
OdeSolver$advance_collect()Advance to each of 'times' in turn and return the state at each: a matrix with 'time' in the first column and one column per state variable, one row per time. The first time must be the current time; the last is landed on exactly. With 'dense = TRUE' the steps are the controller's own and each other time is read off the interpolant of the step that spans it, so the integration costs the same however many rows are asked for. Under '"dopri"' that interpolant has the stepper's own order; under the other two it is cubic Hermite on the step's endpoints, one order short, so prefer '"dopri"' for dense output or 'dense = FALSE', which lands a step on every time.
OdeSolver$advance_collect(times, dense = TRUE)
timesTimes to advance to; the first must be the current time.
dense‘TRUE' to read the requested times off each step’s interpolant, 'FALSE' to land a step on each; see 'advance_collect()'.
OdeSolver$time()Current time.
OdeSolver$time()
OdeSolver$state()Current state.
OdeSolver$state()
OdeSolver$rates()Rates at the current state, as last evaluated (no new evaluation after a step or 'set_state()').
OdeSolver$rates()
OdeSolver$times()Times of every accepted step since the last 'set_state()', starting with the time it was set at.
OdeSolver$times()
OdeSolver$set_state()Reset the state and time. 'y' may have a different length from before. Resets the step size to the control's initial value and the recorded times; evaluates 'rhs' once; clears the effect of an error in a callback.
OdeSolver$set_state(y, time)
yA state vector, of any length.
timeThe time that state is at.
OdeSolver$step_size()The step the controller will try next.
OdeSolver$step_size()
OdeSolver$set_step_size()Set the step the controller tries next, for example to carry a known-good step across a 'set_state()'.
OdeSolver$set_step_size(h)
hA step size.
OdeSolver$counts()What the integration has cost: a named vector with 'n_rhs' (right-hand-side evaluations) and 'n_jac' (Jacobian formations) since the solver was created, 'n_steps' (accepted steps since the last 'set_state()') and 'n_rejections' (attempts rejected and retried smaller, since the solver was created).
OdeSolver$counts()
OdeSolver$clone()The objects of this class are cloneable with this method.
OdeSolver$clone(deep = FALSE)
deepWhether to make a deep clone.
decay <- function(t, y) -y s <- OdeSolver$new(decay, y0 = 1, autonomous = TRUE) s$advance_adaptive(c(0, 1)) s$state() # exp(-1) s$counts()decay <- function(t, y) -y s <- OdeSolver$new(decay, y0 = 1, autonomous = TRUE) s$advance_adaptive(c(0, 1)) s$state() # exp(-1) s$counts()