StackOverFlow Error in SIR model using ModelingTookit

General Usage
, , ,
1

I am brand new to ModelingTookit, and as a first exercise I have defined a simple SIR model and tried to solve it and plot the solution. Unfortunately I’m getting a StackOverFlow Error that I can’t figure out. Here’s my code:

using ModelingToolkit, DifferentialEquations, Plots 

@parameters t β γ N
@variables S(t) I(t) R(t)
D = Differential(t)
eqs = [D(S) ~ -β/N*S*I,
       D(I) ~ β/N*S*I - γ*I,
       D(R) ~ γ*I]

SIR = ODESystem(eqs, name = :SIR)
u0 = [S => 0.99*N,
      I => 0.01*N,
      R => 0.00*N
     ]
p  = [N => 1e3,
      β => 0.20,
      γ => 0.15
      ]
tspan = (0.0,100.0)
prob = ODEProblem(SIR,u0,tspan,p,jac=true)
sol = solve(prob,Tsit5())
plot(sol)

And here’s the stack trace:

StackOverflowError:

Stacktrace:
  [1] recursive_unitless_bottom_eltype(a::Type{Any}) (repeats 65223 times)
    @ RecursiveArrayTools C:\Users\Michael\.julia\packages\RecursiveArrayTools\85a1k\src\utils.jl:91
  [2] recursive_unitless_bottom_eltype(a::Type{Vector{Term{Real}}})
    @ RecursiveArrayTools C:\Users\Michael\.julia\packages\RecursiveArrayTools\85a1k\src\utils.jl:92
  [3] recursive_unitless_bottom_eltype(a::Vector{Term{Real}})
    @ RecursiveArrayTools C:\Users\Michael\.julia\packages\RecursiveArrayTools\85a1k\src\utils.jl:90
  [4] __init(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq C:\Users\Michael\.julia\packages\OrdinaryDiffEq\2ZBfC\src\solve.jl:150
  [5] __init(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq C:\Users\Michael\.julia\packages\OrdinaryDiffEq\2ZBfC\src\solve.jl:67
  [6] __solve(::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq C:\Users\Michael\.julia\packages\OrdinaryDiffEq\2ZBfC\src\solve.jl:4
  [7] __solve(::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5)
    @ OrdinaryDiffEq C:\Users\Michael\.julia\packages\OrdinaryDiffEq\2ZBfC\src\solve.jl:4
  [8] solve_call(_prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:61
  [9] solve_call(_prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:48
 [10] solve_up(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Term{Real}}, p::Vector{Float64}, args::Tsit5; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:82
 [11] solve_up(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Term{Real}}, p::Vector{Float64}, args::Tsit5)
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:75
 [12] solve(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:70
 [13] solve(prob::ODEProblem{Vector{Term{Real}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#200"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0x036c6b2f, 0x3d139a3e, 0xe37bbf78, 0xc69c8058, 0xd46152e3)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTIIPVar#303"), Symbol("##MTKArg#299"), Symbol("##MTKArg#300"), Symbol("##MTKArg#301")), ModelingToolkit.var"#_RGF_ModTag", (0xbf13251e, 0x433c4ae0, 0x9b2a464f, 0x041ba8c5, 0x099b3a9d)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
    @ DiffEqBase C:\Users\Michael\.julia\packages\DiffEqBase\lULzQ\src\solve.jl:68
 [14] top-level scope
    @ In[8]:21
 [15] eval
    @ .\boot.jl:360 [inlined]
 [16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094

Have I defined the model incorrectly? I’m afraid I don’t understand the error…

Also, if I may ask a follow-up question for once the above error is resolved: I’m told that it’s possible to have MTK automatically scale the ODE variables symbolically. For example, in the SIR example I’d like to rewrite the system in terms of the new variables s(t) = S(t)/N, i(t) = I(t)/N, r(t) = R(t)/N, and then solve the system in terms of s, i, r and then plot the solution. How would I implement this?

2

What is typeof(prob.u0) and typeof(prob.p)?

3

typeof(prob.u0) = Vector{Term{Real}} (alias for Array{Term{Real}, 1}).
typeof(prob.p) = Vector{Float64} (alias for Array{Float64, 1})

4

Ohhh, it works now! I changed replaced N with its value in u0:

u0 = [S => 0.99*1e3,
      I => 0.01*1e3,
      R => 0.00*1e3
     ]

So I guess we can’t use variables in the dictionary for u0 then…?

5

N was symbolic so you got a symbolic initial condition which it did not know how to handle. For any numerical solve you need a numerical initial condition.

6

Thanks! Good to know.

Now about my scaling question…how would I go about that? I just want to scale S, I, R by N before solving.

7

Generate a problem with a given N and use remake to change the u0 around.

1 Like
8

Or maybe using substitute before solve is a nice way to do it.

9

So I added s(t) i(t) r(t) to the variables list and then tried

substitute(eqs, Dict([S => s*N, I => i*N, R => r*N]))

but this didn’t seem to change the equations at all:

3-element Vector{Equation}:
 Equation(Differential(S(t)), (((-β) / N) * S(t)) * I(t))
 Equation(Differential(I(t)), (((β / N) * S(t)) * I(t)) - (γ * I(t)))
 Equation(Differential(R(t)), γ * I(t))

I want it to ultimately return

3-element Vector{Equation}:
 Equation(Differential(s(t)), (-β * s(t) * i(t))
 Equation(Differential(i(t)), ((β * s(t)) * i(t)) - (γ * i(t)))
 Equation(Differential(r(t)), γ * i(t))

I also tried substitute([-β/N*S*I], Dict([S => s*N, I => i*N, R => r*N])), but again it didn’t change anything…

10

substitute.(eqs, Dict([S => s*N, I => i*N, R => r*N]))

11

I tried that too but got this error:

ArgumentError: broadcasting over dictionaries and `NamedTuple`s is reserved

Stacktrace:
 [1] broadcastable(#unused#::Dict{Num, Num})
   @ Base.Broadcast .\broadcast.jl:683
 [2] broadcasted(::Function, ::Vector{Equation}, ::Dict{Num, Num})
   @ Base.Broadcast .\broadcast.jl:1313
 [3] top-level scope
   @ In[24]:3
 [4] eval
   @ .\boot.jl:360 [inlined]
 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1094
12

substitute.(eqs, (Dict([S => s*N, I => i*N, R => r*N]),))

2 Likes
13

Thanks! So I now have the following:

eqs = substitute.(eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
eqs = expand_derivatives.(eqs)

s_eq = Symbolics.solve_for(eqs[1],D(s))
i_eq = Symbolics.solve_for(eqs[2],D(i))
r_eq = Symbolics.solve_for(eqs[3],D(r))

eqs_scaled = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]

Which gives me:

\begin{align} \frac{ds(t)}{dt} =& - \beta i\left( t \right) s\left( t \right) \\ \frac{di(t)}{dt} =& - \frac{N \gamma i\left( t \right) - N \beta i\left( t \right) s\left( t \right)}{N} \\ \frac{dr(t)}{dt} =& \gamma i\left( t \right) \end{align}

Now for some reason it doesn’t simplify i(t) by canceling out the N’s. I tried simplify(i_eq) but it just gives back the same thing…is there any way to fully simplify it?

14

Can you make an MWE out of that? I would assume that would further simplify.

15

I managed to figure it out! Here’s a MWE:

using ModelingToolkit, Symbolics

@parameters t β γ N
@variables s(t) i(t) r(t)

D = Differential(t)
eq = D(i) ~ (N*β*s*i - N*γ*i)/(N)     
eq = Symbolics.solve_for(eq,D(i))
Symbolics.simplify(eq)

produces

\begin{equation} \frac{N \beta i\left( t \right) s\left( t \right) - N \gamma i\left( t \right)}{N} \end{equation}

Changing the last line to Symbolics.simplify(eq;expand=true) makes it cancel out the N:
\begin{equation} \beta i\left( t \right) s\left( t \right) - \gamma i\left( t \right) \end{equation}

1 Like
16

@shashi is there a reason why it doesn’t expand that by default?

17

Almost there… here’s what I have so far:

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) 

#1) Define the (unscaled equations)
D = Differential(t)
SIR_eqs = [D(S) ~ -β/N*S*I,
           D(I) ~ β/N*S*I - γ*I,
           D(R) ~ γ*I]

SIR = ODESystem(SIR_eqs, name = :SIR)

#2) Rewrite equations in terms of s = S/N, i = I/N, r = R/N
sir_eqs = substitute.(SIR_eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
sir_eqs = expand_derivatives.(sir_eqs)

s_eq = Symbolics.solve_for(eqs_new[1],D(s))
i_eq = Symbolics.solve_for(eqs_new[2],D(i))
r_eq = Symbolics.solve_for(eqs_new[3],D(r))

#cancel out N 
i_eq = Symbolics.simplify(i_eq;expand=true) 

#Define the new scaled system
sir_eqs = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]
SIR_scaled = ODESystem(eqs_new, name = :SIR_scaled) 

#3) Specify the initial condition and parameter values
u0 = [S => 0.99,
      I => 0.01,
      R => 0.00
     ]

#N omitted from `pars` b.c. it doesn't appear in the scaled system.
pars = [β => 0.20, γ => 0.15]   

#4) Scale the initial condition 
function scale_u0(u0::Vector{Pair{Num, Float64}}, x::Vector{Float64})
    #x = the scaling factors (constants)
    vars = first.(u0)
    vals = last.(u0)
    
    u0_scaled = [Pair(vars[i], x[i]*vals[i]) for i=1:length(vars)]
    return u0_scaled
end    

N = 1e3
u0_scaled = scale_u0(u0,repeat([N],3))

#4) Now create the ODEProblem and solve the scaled system
prob = ODEProblem(SIR_scaled,u0_scaled,(0.0,200.0),pars)
#sol = solve(prob,Tsit5())

But when I try to create the ODEProblem, it gives the following error:

ArgumentError: Term{Real, Nothing}[s(t), i(t), r(t)] are missing from the variable map.

Stacktrace:
  [1] throw_missingvars(vars::Vector{Term{Real, Nothing}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:71
  [2] _varmap_to_vars(varmap::Dict{Num, Float64}, varlist::Vector{Term{Real, Nothing}}; defaults::Dict{Any, Any}, check::Bool, toterm::typeof(Symbolics.diff2term))
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:63
  [3] varmap_to_vars(varmap::Vector{Pair{Num, Float64}}, varlist::Vector{Term{Real, Nothing}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function)
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:38
  [4] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, parammap::Vector{Pair{Num, Float64}}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:461
  [5] process_DEProblem
    @ C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:437 [inlined]
  [6] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:550
  [7] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:550
  [8] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:529
  [9] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N)
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:529
 [10] top-level scope
    @ In[156]:46
 [11] eval
    @ .\boot.jl:360 [inlined]
 [12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094

A quick search brought this up: Missing from the variable map error - Usage - JuliaLang. Not sure if mine is exactly the same issue. I wasn’t able to figure it out…

18

I think you need to specify initial values u0 for lower case [s, i, r] not upper case since you did a substitution. Those are the values mentioned in the error message, the variable map is the list of initial values you provided.

2 Likes
19

Thanks! That was indeed the issue!

20

Alright, here is my completed little toy problem!

using ModelingToolkit, DifferentialEquations, Symbolics

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) x y 

#1) Define the (unscaled equations)
D = Differential(t)
SIR_eqs = [D(S) ~ -β/N*S*I,
           D(I) ~ β/N*S*I - γ*I,
           D(R) ~ γ*I]

SIR = ODESystem(SIR_eqs, name = :SIR)

#2) Rewrite equations in terms of s = S/N, i = I/N, r = R/N
sir_eqs = substitute.(SIR_eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
sir_eqs = expand_derivatives.(sir_eqs)

s_eq = Symbolics.solve_for(sir_eqs[1],D(s))
i_eq = Symbolics.solve_for(sir_eqs[2],D(i))
r_eq = Symbolics.solve_for(sir_eqs[3],D(r))

i_eq = Symbolics.simplify(i_eq;expand=true)   #cancel out the N 


sir_eqs = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]
SIR_scaled = ODESystem(sir_eqs, name = :SIR_scaled)  #Define the new scaled system

#3) Specify the initial condition and parameter values
u0 = [S => 990.0,
      I => 10.0,
      R => 0.00
     ]

pars = [β => 0.20, γ => 0.15]   #N is omitted because it doesn't appear in the scaled system.


function scale_u0(u0::Vector{Pair{Num, Float64}},vars_new::Vector{Num},x::Vector{Float64})
    
    #x = the scaling factors for the variables
    vars = first.(u0)
    vals = last.(u0)
    
    u0_scaled = [Pair(vars_new[i], x[i]*vals[i]) for i=1:length(vars_new)]
    return u0_scaled
end    

N = 1e-3
u0_scaled = scale_u0(u0,[s,i,r],repeat([N],3))

#4) Now solve the ODE and plot it
prob = ODEProblem(SIR_scaled,u0_scaled,(0.0,200.0),pars)
sol = solve(prob,Tsit5(), saveat = 1:1:200)  #works!

Thanks for the help, @ChrisRackauckas! (And also @contradict!).

One follow-up question: Suppose I want to write a function or two that generalizes the above process of scaling the variables to any ODE system whose equations are polynomials in the dependent variables. Ideally I’d like some function that takes in a list of ODE equations, and a vector of variable-scale factor pairs (e.g. [S => 1e3, I => 1e3, R => 1e3], though the scaling factors don’t need to be the same for each). Then I’d want the function to return a new ODE system written in terms of the scaled variables. Is this possible the way that ModelingToolkit is written? I something like this could ever be added to the MKT package, that would really be awesome.