The odeintr is package for integrating differential equations in R. The integration engine is the Boost odeint package.
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.
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()
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.
Pull requests are welcome.