implicit: Compile Stiff ODE system solver

In odeintr: C++ ODE Solvers Compiled on-Demand

Description

Generates a stiff integrator using Rcpp

Usage

compile_implicit(name, sys, pars = NULL, const = FALSE, sys_dim = -1L,
  jacobian = JacobianCpp(sys, sys_dim), atol = 1e-06, rtol = 1e-06,
  globals = "", headers = "", footers = "", compile = TRUE,
  env = new.env(), ...)

JacobianCpp(sys, sys_dim = -1L)

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
sys_dim	length of the state vector
jacobian	C++ expressions computing the Jacobian
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
env	install functions into this environment
...	passed to sourceCpp

Details

See compile_sys for details. This function behaves identically except that you cannot specify a custom observer and you must provide a Jacobian C++ function body. By default, the Jacobian will be symbollically computed from the system function using the JacobianCpp function. This uses D to compute the symbolic derivatives. The integration methods is the 4th-order Rosenbrock implicit solver. Step size is adaptive and output is interpolated.

If you provide a hand-written Jacobian, it must update the matrix J and the vector dfdt. It is perhaps easiest to use JacobianCpp to see the required format.

Author(s)

Timothy H. Keitt

See Also

set_optimization, compile_sys

Examples

## Not run: 
# 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
cat(JacobianCpp(Lorenz.sys))
compile_implicit("lorenz", Lorenz.sys, pars, TRUE)
x = lorenz(rep(1, 3), 100, 0.001)
plot(x[, c(2, 4)], type = 'l', col = "steelblue")
Sys.sleep(0.5)
# Stiff example from odeint docs
Stiff.sys = '
 dxdt[0] = -101.0 * x[0] - 100.0 * x[1];
 dxdt[1] = x[0];
' # Stiff.sys
cat(JacobianCpp(Stiff.sys))
compile_implicit("stiff", Stiff.sys)
x = stiff(c(2, 1), 7, 0.001)
plot(x[, 1:2], type = "l", lwd = 2, col = "steelblue")
lines(x[, c(1, 3)], lwd = 2, col = "darkred")
# Robertson chemical kinetics problem
Robertson = '
dxdt[0] = -alpha * x[0] + beta * x[1] * x[2];
dxdt[1] = alpha * x[0] - beta * x[1] * x[2] - gamma * x[1] * x[1];
dxdt[2] = gamma * x[1] * x[1];
' # Robertson
pars = c(alpha = 0.04, beta = 1e4, gamma = 3e7)
init.cond = c(1, 0, 0)
cat(JacobianCpp(Robertson))
compile_implicit("robertson", Robertson, pars, TRUE)
at = 10 ^ seq(-5, 5, len = 400)
x = robertson_at(init.cond, at)
par(mfrow = c(3, 1), mar = rep(0.5, 4), oma = rep(5, 4), xpd = NA)
plot(x[, 1:2], type = "l", lwd = 3,
     col = "steelblue", log = "x", axes = F, xlab = NA)
axis(2); box()
plot(x[, c(1, 3)], type = "l", lwd = 3,
     col = "steelblue", log = "x", axes = F, xlab = NA)
axis(4); box()
plot(x[, c(1, 4)], type = "l", lwd = 3,
     col = "steelblue", log = "x", axes = F)
axis(2); axis(1); box()

## End(Not run)

odeintr documentation built on May 2, 2019, 2:08 a.m.