targeted::ode_solve: Solving Ordinary Differential Equations

Klaus Kähler Holst

2024-11-25

Introduction

Mathematical and statistical software often relies on sequential computations. ordinary differential equations where numerical approximations are based on looping over the evolving time. When using high-level languages such as R calculations can be very slow unless the algorithms can be vectorized

There are various excellent (O)DE solvers (R: deSolve) Here I will illustrate the above techniques using the targeted R-package based on the target C++ library.

The ODE is specified using the specify_ode function

args(targeted::specify_ode)
#> function (code, fname = NULL, pname = c("dy", "x", "y", "p")) 
#> NULL

The differential equations are here specified as a string containing the C++ code defining the differential equation via the code argument. The variable names are defined through the pname argument which defaults to

dy

Vector with derivatives, i.e. the rhs of the ODE, $y'(t)$ (the result).

x

Vector, with the first element being the time, and the following elements additional exogenous input variables, $x(t) = \{t, x_{1}(t), \ldots, x_{k}(t)\}$

y

Dependent variable, $y(t) = \{y_{1}(t),\ldots,y_{l}(t)\}$

p

Parameter vector $$y'(t) = f_{p}(x(t), y(t))$$

All variables are treated as armadillo vectors/matrices, arma::mat.

As an example, we can specify the simple differential equation $$y'(t) = y(t)-1$$

dy <- targeted::specify_ode("dy = y - 1;")

This compiles the function and stores the pointer in the variable dy.

To solve the ODE we must then use the function solve_ode

args(targeted::solve_ode)
#> function (ode_ptr, input, init, par = 0) 
#> NULL

The first argument is the external pointer, the second argument input is the input matrix ($x(t)$ above), and the init argument is the vector of initial boundary conditions $y(0)$. The argument par is the vector of parameters defining the ODE ($p$).

Examples

In this example the input variable does not depend on any exogenous variables so we only need to supply the time points, and the defined ODE does not depend on any parameters. To approximate the solution with initial condition $y(0)=0$, we therefore run the following code

t <- seq(0, 10, length.out=1e4)
y <- targeted::solve_ode(dy, t, init=0)
plot(t, y, type='l', lwd=3)

As a more interesting example consider the Lorenz Equations $$\frac{dx_{t}}{dt} = \sigma(y_{t}-x_{t})$$$$\frac{dy_{t}}{dt} = x_{t}(\rho-z_{t})-y_{t}$$$$\frac{dz_{t}}{dt} = x_{t}y_{t}-\beta z_{t}$$

we may define them as

library(targeted)
ode <- 'dy(0) = p(0)*(y(1)-y(0));
        dy(1) = y(0)*(p(1)-y(2));
        dy(2) = y(0)*y(1)-p(2)*y(2);'
f <- specify_ode(ode)

With the choice of parameters given by $\sigma=10, \rho=28, \beta=8/3$ and initial conditions $(x_0,y_0,z_0)=(1,1,1)$, we can calculate the solution

tt <- seq(0, 100, length.out=2e4)
y <- solve_ode(f, input=tt, init=c(1, 1, 1), par=c(10, 28, 8/3))
head(y)
#>          [,1]     [,2]      [,3]
#> [1,] 1.000000 1.000000 1.0000000
#> [2,] 1.003322 1.135177 0.9920639
#> [3,] 1.013094 1.271279 0.9849468
#> [4,] 1.029068 1.409156 0.9786954
#> [5,] 1.051048 1.549630 0.9733721
#> [6,] 1.078888 1.693496 0.9690541

colnames(y) <- c("x","y","z")
scatterplot3d::scatterplot3d(y, cex.symbols=0.1, type='b',
                             color=viridisLite::viridis(nrow(y)))

To illustrate the use of exogenous inputs, consider the following simulated data

n <- 1e4
tt <- seq(0, 10, length.out=n)  # Time
xx <- rep(0, n)
xx[(n / 3):(2 * n / 3)] <- 1 # Exogenous input, x(t)
input <- cbind(tt, xx)

and the following ODE

$$y'(t) = \beta_{0} + \beta_{1}y(t) + \beta_{2}y(t)x(t) + \beta_{3}x(t)\cdot t$$

mod <- 'double t = x(0);
        dy = p(0) + p(1)*y + p(2)*x(1)*y + p(3)*x(1)*t;'
dy <- specify_ode(mod)

With $y(0)=100$ and $\beta_0=0, \beta_{1}=0.4, \beta_{2}=-0.5, \beta_{3}=-5$ we obtain the following solution

yy <- solve_ode(dy, input=input, init=100, c(0, .4, -.5, -5))
plot(tt, yy, type='l', lwd=3, xlab='time', ylab='y')

SessionInfo

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] targeted_0.6 lava_1.8.0  
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-1           future.apply_1.11.3    jsonlite_1.8.9        
#>  [4] compiler_4.4.2         Rcpp_1.0.13-1          parallel_4.4.2        
#>  [7] jquerylib_0.1.4        globals_0.16.3         splines_4.4.2         
#> [10] systemfonts_1.1.0      textshaping_0.4.0      yaml_2.3.10           
#> [13] fastmap_1.2.0          lattice_0.22-6         R6_2.5.1              
#> [16] knitr_1.49             future_1.34.0          nloptr_2.1.1          
#> [19] desc_1.4.3             bslib_0.8.0            rlang_1.1.4           
#> [22] cachem_1.1.0           xfun_0.49              fs_1.6.5              
#> [25] sass_0.4.9             viridisLite_0.4.2      cli_3.6.3             
#> [28] pkgdown_2.1.1          digest_0.6.37          grid_4.4.2            
#> [31] lifecycle_1.0.4        RcppArmadillo_14.2.0-1 scatterplot3d_0.3-44  
#> [34] evaluate_1.0.1         pracma_2.4.4           data.table_1.16.2     
#> [37] numDeriv_2016.8-1.1    listenv_0.9.1          codetools_0.2-20      
#> [40] ragg_1.3.3             survival_3.7-0         optimx_2023-10.21     
#> [43] parallelly_1.39.0      rmarkdown_2.29         tools_4.4.2           
#> [46] htmltools_0.5.8.1