compile_sys: Compile ODE system
In odeintr: C++ ODE Solvers Compiled on-Demand

Description Usage Arguments Details Value Note Author(s) See Also Examples

View source: R/odeintr.R

Generates an integrator using Rcpp

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(), ...)

`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`

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
`bsd`	`bulirsch_stoer_dense_out`	Interpolating adaptive

These steppers are described at here.

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
`name_set_params`	set parameter values	parameters	void

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

The c++11 plugin is enabled.

Timothy H. Keitt

set_optimization, integrate_sys

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

## End(Not run)