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
library(odeintr) dxdt = function(x, t) x * (1 - x) system.time({x = integrate_sys(dxdt, 0.001, 15, 0.01)}) 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)}) 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)}) 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)}) 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)}) 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)) 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")
The example below is not working with the latest version of odeint that comes with the BH package. I've submitted and issue on github.
# 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()
It is important to understand that the only thing that odeintr
does is to generate C++ code and compile it using Rcpp. That means you can use the generated code as a starting point for your project and modify as you wish.
the_code = compile_sys("logitest", "dxdt = x * (1 - x)", compile = FALSE) cat(the_code)
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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.