Numerical ODE Solver

odeintr by Timothy H. Keitt is an R package for integrating differential equations with the Boost odeint library. It allows specifying the model in C++ code and compiling it on the fly with Rcpp.

Here we provide an interface to the odeintr package, but

  1. Use an intuitive way to specify the ODE model with symengine objects;
  2. Use symengine’s code generation functionality to generate the C++ source and compile with odeintr.

However, the interface is not stable and may subject to change in future. This mainly serves an example of how to extend symengine and you are welcome to check out the source code.

library(symengine)
## SymEngine Version: 0.11.2
##  _____           _____         _         
## |   __|_ _ _____|   __|___ ___|_|___ ___ 
## |__   | | |     |   __|   | . | |   | -_|
## |_____|_  |_|_|_|_____|_|_|_  |_|_|_|___|
##       |___|               |___|

Define ordinary differential equations with symengine

A ordinary differential equation could be constructed with dxdt(x) == rhs where x and rhs will be converted to SymEngine’s ‘Basic’ S4 object. This works by defining a S4 method of == for the return type of dxdt().

For example

x <- Symbol("x")
a <- 3
eq <- dxdt(x) == a/(x + 1)
eq
## Ordinary differential equation:
## d(x)/dt == 3.0/(1.0 + x)

Define ODE system

ODESystem will take a list of ordinary differential equations, generate C++ source code and compile on the fly with Rcpp. The following is the Lorenz system.

sigma <- 10
rho <- 28
beta <- 8/3
use_vars(x, y, z)
## Initializing 'x', 'y', 'z'
lorenz_sys <- list(
    dxdt(x) == sigma * (y - x),
    dxdt(y) == (rho - z) * x - y,
    dxdt(z) == - beta * z + x * y
)
lorenz_sys <- ODESystem(lorenz_sys, method = "rk5_i")

The method argument is passed to odeintr::compile_sys.

Get results

A S4 method of predict is defined to run the model with given initial conditions, duration and step_size.

res <- predict(lorenz_sys, init = c(x=1, y=1, z=1),
               duration = 100, step_size = 0.001, start = 0)
head(res)
##    Time        x        y         z
## 1 0.000 1.000000 1.000000 1.0000000
## 2 0.001 1.000130 1.025989 0.9983486
## 3 0.002 1.000516 1.051961 0.9967278
## 4 0.003 1.001158 1.077922 0.9951378
## 5 0.004 1.002051 1.103879 0.9935788
## 6 0.005 1.003193 1.129840 0.9920511
plot(res[, c(2, 4)], type = 'l', col = "steelblue", main = "Lorenz Attractor")

Van der Pol Oscillator

Example of Van der Pol oscillator from odeintr package.

use_vars(x, y)
## Initializing 'x', 'y'
vdp_sys <- ODESystem(
    dxdt(x) == y,
    dxdt(y) == 2 * (1 - x * x) * y - x,
    method = "bsd" # Bulirsch-Stoer
)
res <- predict(vdp_sys, init = rep(1e-4, 2), duration = 100, step_size = 0.01)
oldpar <- 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(res[, c(3, 1)]); axis(3); axis(4)
make.plot(res[, c(1, 2)], "Time", "X1"); axis(1); axis(2)
make.plot(res[, c(3, 2)], "X2"); axis(1); axis(4)
title(main = "Van der Pol Oscillator", outer = TRUE)

par(oldpar)