Travis-CI Build Status CRAN Version CRAN Downloads

The odeintr is package for integrating differential equations in R. The integration engine is the Boost odeint package.

Features

  1. Simple specification of the ODE system
  2. Named, dynamic run-time setable system parameters
  3. Intelligent defaults, easily overridden, used throughout
  4. A wide range of integration methods available for compiled system (see stepper types)
  5. Fully automated compilation of ODE system specified in C++
  6. Simple openmp vectorization of large systems (Broken in latest Boost release)
  7. Results returned as a simple data frame ready for analysis and plotting
  8. Ability to specify a custom observer in R that can return arbitrary data
  9. Three options for calling the observer: at regular intervals, after each update step or at specified times
  10. Ability to alter system state and restart simulations where you left off
  11. Can compile an implicit solver with symbolic evaluation of the Jacobian
  12. You can easily save and edit the generated C++ code

Installation

install.packages(odeintr)                   # released
devtools::install_github("thk686/odeintr")  # development

There have been some problems compiling on Windows owing to an incompatibility between the BH package and RTools toolchain. I have removed the dependency on the BH package. If you have trouble with Windows, please try the development version.

Examples

library(odeintr)
dxdt = function(x, t) x * (1 - x)
system.time({x = integrate_sys(dxdt, 0.001, 15, 0.01)})
##    user  system elapsed 
##   0.093   0.016   0.109
plot(x, type = "l", lwd = 3, col = "steelblue", main = "Logistic Growth")

compile_sys("logistic", "dxdt = x * (1 - x)")
system.time({x = logistic(0.001, 15, 0.01)})
##    user  system elapsed 
##   0.000   0.000   0.001
plot(x, type = "l", lwd = "3", col = "steelblue", main = "Logistic Growth")
points(logistic_at(0.001, sort(runif(10, 0, 15)), 0.01), col = "darkred")

dxdt = function(x, t) c(x[1] - x[1] * x[2], x[1] * x[2] - x[2])
obs = function(x, t) c(Prey = x[1], Predator = x[2], Ratio = x[1] / x[2])
system.time({x = integrate_sys(dxdt, rep(2, 2), 20, 0.01, observer = obs)})
##    user  system elapsed 
##   0.295   0.042   0.473
plot(x[, c(2, 3)], type = "l", lwd = 2, col = "steelblue", main = "Lotka-Volterra Phase Plot")

plot(x[, c(1, 4)], type = "l", lwd = 2, col = "steelblue", main = "Prey-Predator Ratio")

# C++ code
Lorenz.sys = '
  dxdt[0] = 10.0 * (x[1] - x[0]);
  dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
  dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
  ' # Lorenz.sys
compile_sys("lorenz", Lorenz.sys)
system.time({x = lorenz(rep(1, 3), 100, 0.001)})
##    user  system elapsed 
##   0.009   0.012   0.033
plot(x[, c(2, 4)], type = 'l', col = "steelblue", main = "Lorenz Attractor")

VdP.sys = '
dxdt[0] = x[1];
dxdt[1] = 2.0 * (1 - x[0] * x[0]) * x[1] - x[0];
' # Vdp.sys
compile_sys("vanderpol", VdP.sys, method = "bsd") # Bulirsch-Stoer
system.time({x = vanderpol(rep(1e-4, 2), 100, 0.01)})
##    user  system elapsed 
##   0.004   0.000   0.003
par(mfrow = c(2, 2), mar = rep(0.5, 4), oma = rep(5, 4), xpd = NA)
make.plot = function(xy, xlab = NA, ylab = NA)
  plot(xy, col = "steelblue", lwd = 2, type = "l",
       axes = FALSE, xlab = xlab, ylab = ylab)
plot.new()
make.plot(x[, c(3, 1)]); axis(3); axis(4)
make.plot(x[, c(1, 2)], "Time", "X1"); axis(1); axis(2)
make.plot(x[, c(3, 2)], "X2"); axis(1); axis(4)
title(main = "Van der Pol Oscillator", outer = TRUE)

# Use a dynamic parameter
VdP.sys = '
dxdt[0] = x[1];
dxdt[1] = mu * (1 - x[0] * x[0]) * x[1] - x[0];
' # Vdp.sys
compile_sys("vpol2", VdP.sys, "mu", method = "bsd")
par(mfrow = c(2, 2), mar = rep(1, 4), oma = rep(3, 4), xpd = NA)
for (mu in seq(0.5, 2, len = 4))
{
  vpol2_set_params(mu = mu)
  x = vpol2(rep(1e-4, 2), 100, 0.01)
  make.plot(x[, 2:3]); box()
  title(paste("mu =", round(mu, 2)))
}
title("Van der Pol Oscillator Parameter Sweep", outer = TRUE)
title(xlab = "X1", ylab = "X2", line = 0, outer = TRUE)

# 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))
## J(0, 0) = -101;
## J(0, 1) = -100;
## J(1, 0) = 1;
## J(1, 1) = 0;
## dfdt[0] = 0.0;
## dfdt[1] = 0.0;
compile_implicit("stiff", Stiff.sys)
x = stiff(c(2, 1), 5, 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))
## J(0, 0) = -alpha;
## J(0, 1) = beta * x[2];
## J(0, 2) = beta * x[1];
## J(1, 0) = alpha;
## J(1, 1) = -(beta * x[2] + (gamma * x[1] + gamma * x[1]));
## J(1, 2) = -(beta * x[1]);
## J(2, 0) = 0;
## J(2, 1) = gamma * x[1] + gamma * x[1];
## J(2, 2) = 0;
## dfdt[0] = 0.0;
## dfdt[1] = 0.0;
## dfdt[2] = 0.0;
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()

Performance

Because ODEINT is a header-only library, the entire integration path is exposed to the compiler. That means your system functions can be inlined with the integration code, loops unrolled, etc. It will help if you enable optimization in your compiler. Use “-O3” with gcc. See the R documentation on the user Makevars file. (Odeintr now provides a convenient function to set the compiler optimization level.)

The Lorenz and Van der Pol examples above show about 10 million observer calls per second.

To Do

  1. Add additional integration methods from odeint
  2. Extend customized observer to compiled code
  3. Allow user to set error tolerances in compiled code
  4. Allow user to set error tolerances for system defined in R
  5. Expose implicit solver methods
  6. Compute Jacobian symbolically
  7. Convenient dynamic parameter settings
  8. Install emitted function in a new environment
  9. Add more control of openmp threading
  10. Allow custom state type

Pull requests are welcome.