Solving Ordinary Differential Equations (ODE) in R with diffeqr

1D Linear ODEs

Let’s solve the linear ODE u'=1.01u. First setup the package:

diffeqr::diffeq_setup()

Define the derivative function f(u,p,t).

f <- function(u,p,t) {
  return(1.01*u)
}

Then we give it an initial condition and a time span to solve over:

u0 = 1/2
tspan <- list(0.0,1.0)

With those pieces we call diffeqr::ode.solve to solve the ODE:

sol = diffeqr::ode.solve(f,u0,tspan)

This gives back a solution object for which sol$t are the time points and sol$u are the values. We can check it by plotting the solution:

plot(sol$t,sol$u,"l")

linear_ode

Systems of ODEs

Now let’s solve the Lorenz equations. In this case, our initial condition is a vector and our derivative functions takes in the vector to return a vector (note: arbitrary dimensional arrays are allowed). We would define this as:

f <- function(u,p,t) {
  du1 = p[1]*(u[2]-u[1])
  du2 = u[1]*(p[2]-u[3]) - u[2]
  du3 = u[1]*u[2] - p[3]*u[3]
  return(c(du1,du2,du3))
}

Here we utilized the parameter array p. Thus we use diffeqr::ode.solve like before, but also pass in parameters this time:

u0 = c(1.0,0.0,0.0)
tspan <- list(0.0,100.0)
p = c(10.0,28.0,8/3)
sol = diffeqr::ode.solve(f,u0,tspan,p=p)

The returned solution is like before. It is convenient to turn it into a data.frame:

udf = as.data.frame(sol$u)

Now we can use matplot to plot the timeseries together:

matplot(sol$t,udf,"l",col=1:3)

Now we can use the Plotly package to draw a phase plot:

plotly::plot_ly(udf, x = ~V1, y = ~V2, z = ~V3, type = 'scatter3d', mode = 'lines')

timeseries

Plotly is much prettier!

Option Handling

If we want to have a more accurate solution, we can send abstol and reltol. Defaults are 1e-6 and 1e-3 respectively. Generally you can think of the digits of accuracy as related to 1 plus the exponent of the relative tolerance, so the default is two digits of accuracy. Absolute tolernace is the accuracy near 0.

In addition, we may want to choose to save at more time points. We do this by giving an array of values to save at as saveat. Together, this looks like:

abstol = 1e-8
reltol = 1e-8
saveat = 0:10000/100
sol = diffeqr::ode.solve(f,u0,tspan,p=p,abstol=abstol,reltol=reltol,saveat=saveat)
udf = as.data.frame(sol$u)
plotly::plot_ly(udf, x = ~V1, y = ~V2, z = ~V3, type = 'scatter3d', mode = 'lines')

precise_solution

We can also choose to use a different algorithm. The choice is done using a string that matches the Julia syntax. See the ODE tutorial for details. The list of choices for ODEs can be found at the ODE Solvers page. For example, let’s use a 9th order method due to Verner:

sol = diffeqr::ode.solve(f,u0,tspan,alg="Vern9()",p=p,abstol=abstol,reltol=reltol,saveat=saveat)

Note that each algorithm choice will cause a JIT compilation.

Performance Enhancements

One way to enhance the performance of your code is to define the function in Julia so that way it is JIT compiled. diffeqr is built using the JuliaCall package, and so you can utilize this same interface to define a function directly in Julia:

f <- JuliaCall::julia_eval("
function f(du,u,p,t)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end")

We can then use this in our ODE function by telling it to use the Julia-defined function called f:

u0 = c(1.0,0.0,0.0)
tspan <- list(0.0,100.0)
sol = diffeqr::ode.solve('f',u0,tspan)

This will help a lot if you are solving difficult equations (ex. large PDEs) or repeat solving (ex. parameter estimation).