Nicer way to plug an R function into the RHS of an ODEProblem?

Specific Domains Modelling & Simulations
1

Hi All,

I’m trying to plug an R function into an ODEProblem in OrdinaryDiffEq. The below code works, but is there a simple way to make the code neater e.g. define ODEProblem in terms of an RObject?

using OrdinaryDiffEq
using RCall
using Plots

R"""
sir_ode_r <- function(u,p,t){
    S <- u[1]
    I <- u[2]
    R <- u[3]
    N <- S+I+R
    beta <- p[1]
    cee <- p[2]
    gamma <- p[3]
    dS <- -beta*cee*I/N*S
    dI <- beta*cee*I/N*S - gamma*I
    dR <- gamma*I
    return(c(dS,dI,dR))
}
"""

function sir_ode_jl(u,p,t)
    robj = rcall(:sir_ode_r, u, p, t)
    return convert(Array,robj)
end

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem(sir_ode_jl, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)
1 Like
2

For this kind of function, you can use ModelingToolkit to translate it to a symbolic form and generate the Julia function. That would also fix the performance. That’s how the diffeqr/diffeqpy JIT compiler works.

Otherwise, I think this is how I would expect it to look.

3

Cool - can you show a simple example (modeltoolkitize?)

Would this also work with my PythonCall example?

using OrdinaryDiffEq
using PythonCall
using Plots

@pyexec """
def sir_ode_py(u,p,t):
    S = u[0]
    I = u[1]
    R = u[2]
    N = S+I+R
    beta = p[0]
    c = p[1]
    gamma = p[2]
    dS = -beta*c*I/N*S
    dI = beta*c*I/N*S - gamma*I
    dR = gamma*I
    return [dS,dI,dR]
""" => sir_ode_py

sir_ode_jl(u,p,t) = pyconvert(Array{Float64}, sir_ode_py(u, p, t))

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem{false}(sir_ode_jl, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)
4

It should. Try just calling modelingtoolkitize(prob). If you make that convert into an Array instead of forcing Array{Float64} you should be good.

5

For Python, this code works (converting to Array rather than Array{Float64}

using ModelingToolkit
@named sys = modelingtoolkitize(prob)
prob_mtk = ODEProblem(sys, u0, tspan, p)
sol_mtk = solve(prob_mtk, Tsit5(), dt = δt)
plot(sol_mtk)

However, modelingtoolkitize breaks for the R version: ERROR: MethodError: no method matching sexpclass(::Num)

(yes, I know I’m not using MTK9 :slight_smile: )

6

PS. Is there any way to pass the names of states and parameters when using modelingtoolkitize so that the resulting equations are easier to read?

7

No, but it wouldn’t take more than 15 minutes to add. Open an issue so hopefully I don’t forget. At a conference right now but it’ll end up in the email list that way.