Specify names of states and parameters in ODEProblem?

Specific Domains Modelling & Simulations
1

Here’s a simple problem I have where I couldn’t find any documentation; I’d like to specify the parameter names of the states and parameters to SciML’s ODEFunction e.g. using Dict(:u => [:S, :I, :R], :p => [:β, :c, :γ]). Here’s a complete example that doesn’t specify the names. I know that this is done automatically using ModelingToolkit, but I want to use vanilla Julia. Any suggestions?

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end
prob_ode = ODEProblem(sir_ode!, u0, tspan, p)
sol_ode = solve(prob_ode, Tsit5(), dt = δt)
2

This is all defined in Home · SymbolicIndexingInterface.jl

3

Thanks Chris! Just what I was looking for…

4

One more thing; while this gives me a way to index the solution (e.g. sol[:S]), the names don’t show up automatically when I plot the solution (U just get u[1] etc.). Is there a way to show the defined symbols when plotting?

using OrdinaryDiffEq
using SymbolicIndexingInterface
using ModelingToolkit
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end
sys = SymbolCache([:S,:I,:R],[:β,:c,:γ],:t)
fun = ODEFunction(sir_ode!; sys=sys)
prob = ODEProblem(fun, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)
5

The plots should already use that? If not, then it needs some work in the plot recipe. It just asks the sys for the names to generate the labels. @cryptic.ax is SymbolCache missing a dispatch on nameof or something?

6

The plot recipe definitely doesn’t use that; here’s the full code for debugging purposes:

using OrdinaryDiffEq
using SymbolicIndexingInterface
using ModelingToolkit
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end
sys = SymbolCache([:S,:I,:R],[:β,:c,:γ],:t)
fun = ODEFunction(sir_ode!; sys=sys)

δ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(fun, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)

The symbols are also lost during modelingtoolkitize too:

sys_mtk = modelingtoolkitize(prob)

gives:

julia> sys_mtk.eqs
Equation[Differential(t)(x₁(t)) ~ (-α₁*α₂*x₁(t)*x₂(t)) / (x₁(t) + x₂(t) + x₃(t)), Differential(t)(x₂(t)) ~ (α₁*α₂*x₁(t)*x₂(t)) / (x₁(t) + x₂(t) + x₃(t)) - α₃*x₂(t), Differential(t)(x₃(t)) ~ α₃*x₂(t)][1:3]

Model ##MTKizedODE#710 with 3 equations
States (3):
  x₁(t) [defaults to 990.0]
  x₂(t) [defaults to 10.0]
⋮
Parameters (3):
  α₁ [defaults to 0.05]
  α₂ [defaults to 10.0]
⋮
7

I’ll look into this, the plot recipes should respect the symbol names.

The symbols are also lost during modelingtoolkitize too:

That’s not implemented yet, but it is possible.