compile_sys: Compile ODE system

Description

Generates an integrator using Rcpp

Usage

compile_sys(name, sys, pars = NULL, const = FALSE, method = "rk5_i",
  sys_dim = -1L, atol = 1e-06, rtol = 1e-06, globals = "",
  headers = "", footers = "", compile = TRUE, observer = NULL,
  env = new.env(), ...)

Arguments

name

the name of the generated integration function

sys

a string containing C++ expressions

pars

a named vector of numbers or a vector of names or number of parameters

const

declare parameters const if true

method

a method string (see Details)

sys_dim

length of the state vector

atol

absolute tolerance if using adaptive step size

rtol

relative tolerance if using adaptive step size

globals

a string with global C++ declarations

headers

code to appear before the odeintr namespace

footers

code to appear after the odeintr namespace

compile

if false, just return the code

observer

an optional R function to record output

env

install functions into this environment

...

passed to sourceCpp

Value

The C++ code invisibly. The following functions are generated:

Function	Use	Arguments
Return	`name`	regular observer calls
`init, duration, step_size = 1.0, start = 0.0`	data frame	`name_adap`
adaptive observer calls	`init, duration, step_size = 1.0, start = 0.0`	data frame
`name_at`	specified observer calls	`init, times, step_size = 1.0, start = 0.0`
data frame	`name_continue_at`	specified observer calls starting from previous final state
`times, step_size = 1.0`	data frame	`name_no_record`
no observer calls	`init, duration, step_size = 1.0, start = 0.0`	vector (final state)
`name_reset_observer`	clear observed record	void
void	`name_get_state`	get current state
void	vector	`name_set_state`
set current state	`new_state`	void
`name_get_output`	fetch observed record	void
data frame	`name_get_params`	get parameter values
void	a list	Function

Arguments are:

init

vector of initial conditions

duration

end at start + duration

step_size

the integration step size; variable for adaptive methods

start

the starting time (always equal 0.0 for name_at)

time

vector of times as which to call the observer

new_state

vector of state values

Details

C++ code is generated and compiled with sourceCpp. The returned function will integrate the system starting from a provided initial condition and initial time to a specified final time. An attempt is made to get the length of the state vector from the system definition. If this fails, the code will likely crash your R session. It is safer to specify sys_dim directly. The C++ expressions must index a state array of length sys_dim. The state array is x and the derivatives are dxdt. The first state value is x[0] and the first derivative is dxdt[0]. In the case you use bare dxdt and x, an attempt will be made to append [0] to these variables. This can fail, so do not rely on it. This will also fail if you set sys_dim to a positive value. The globals string can be arbitrary C++ code. You can set global named parameter values here. Note that these will be defined within the odeintr namespace. If you supply the pars argument, these parameters will be compiled into the code. There are three options: 1) if pars is a single number, then you can access a vector of parameters named pars of the specified length; 2) if pars is a character vectors, then a parameter will be defined for each; and 3) if the character vector is named, then the names will be used for the parameter names and the associated values will be used as defaults. If you specify const = TRUE, these named parameters will be declared const. Otherwise parameter getter/setter functions will be defined. If observer is an R function, then this function will be used to record the output of the integration. It is called with signature obsever(x, t). Its return value will be coerced to a list. Observer getter/setter functions will be emitted as well (name_g(s)et_observer). You can also get and set an output processing function (name_g(s)et_output_processor). It will be passed a 2-element list. The first element is a vector of time points and the 2nd element is a list of lists, one list per time point. The default processor converts this to a data frame. You can insert arbitrary code outside the odeintr name space using headers and footers. This code can be anything compatible with Rcpp. You could for example define exported Rcpp functions that set simulation paramters. headers is inserted right after the Rcpp and ODEINT includes. footers is inserted at the end of the code. The following methods can be used:

Code	Stepper
Type	`euler`
`euler`	Interpolating
`rk4`	`runge_kutta4`
Regular	`rk54`
`runge_kutta_cash_karp54`	Regular
`rk54_a`	`runge_kutta_cash_karp54`
Adaptive	`rk5`
`runge_kutta_dopri5`	Regular
`rk5_a`	`runge_kutta_dopri5`
Adaptive	`rk5_i`
`runge_kutta_dopri5`	Interpolating adaptive
`rk78`	`runge_kutta_fehlberg78`
Regular	`rk78_a`
`runge_kutta_fehlberg78`	Adaptive
`abN`	`adams_bashforth`
Order N multistep	`abmN`
`adams_bashforth_moulton`	Order N multistep
`bs`	`bulirsch_stoer`
Adaptive	Code

These steppers are described at https://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.stepper_overview.

Examples

Run this code

## Not run: ------------------------------------
# # Logistic growth
# compile_sys("logistic", "dxdt = x * (1 - x)")
# plot(logistic(0.001, 15, 0.1), type = "l", lwd = 2, col = "steelblue")
# Sys.sleep(0.5)
# 
# # Lotka-Volterra predator-prey equations
# pars = c(alpha = 1, beta = 1, gamma = 1, delta = 1)
# LV.sys = '
#   dxdt[0] = alpha * x[0] - beta * x[0] * x[1];
#   dxdt[1] = gamma * x[0] * x[1] - delta * x[1];
# ' # LV.sys
# compile_sys("preypred", LV.sys, pars, TRUE)
# x = preypred(rep(2, 2), 100, 0.01)
# plot(x[, 2:3], type = "l", lwd = 2,
#      xlab = "Prey", ylab = "Predator",
#      col = "steelblue")
# Sys.sleep(0.5)
# 
# # Lorenz model from odeint examples
# pars = c(sigma = 10, R = 28, b = 8 / 3)
# Lorenz.sys = '
#   dxdt[0] = sigma * (x[1] - x[0]);
#   dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
#   dxdt[2] = -b * x[2] + x[0] * x[1];
# ' # Lorenz.sys
# compile_sys("lorenz", Lorenz.sys, pars, TRUE)
# x = lorenz(rep(1, 3), 100, 0.001)
# plot(x[, c(2, 4)], type = 'l', col = "steelblue")
## ---------------------------------------------