odeintr (version 1.7.0)

compile_implicit: Compile Stiff ODE system solver

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.

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