# StackOverFlow Error in SIR model using ModelingTookit

**URL:** https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604
**Category:** General Usage
**Tags:** error, ode, stackoverflow, modelingtoolkit
**Created:** [July 13, 2021, 9:22pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604 "2021-07-13T21:22:26Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 13, 2021, 9:22pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/1 "2021-07-13T21:22:26Z")

</div>

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:

```julia
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:

```julia
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?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 13, 2021, 9:30pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/2 "2021-07-13T21:30:17Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 13, 2021, 9:42pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/3 "2021-07-13T21:42:19Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 13, 2021, 9:45pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/4 "2021-07-13T21:45:04Z")

</div>

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

```julia
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…?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 13, 2021, 9:50pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/5 "2021-07-13T21:50:49Z")

</div>

> [@Michael\_Barmann](#):
>
> Vector{Term{Real}}

`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.

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 13, 2021, 9:54pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/6 "2021-07-13T21:54:51Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 13, 2021, 9:55pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/7 "2021-07-13T21:55:23Z")

</div>

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

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 13, 2021, 9:56pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/8 "2021-07-13T21:56:09Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 14, 2021, 10:39pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/9 "2021-07-14T22:39:00Z")

</div>

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:

```julia
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

```julia
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…

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 15, 2021, 2:52am UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/10 "2021-07-15T02:52:14Z")

</div>

> [@Michael\_Barmann](#):
>
> substitute(eqs, Dict([S =\> s\*N, I =\> i\*N, R =\> r\*N]))

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 15, 2021, 3:13pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/11 "2021-07-15T15:13:30Z")

</div>

> [@ChrisRackauckas](#):
>
> substitute.(eqs, Dict([S =\> s\*N, I =\> i\*N, R =\> r\*N]))

I tried that too but got this error:

```julia
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

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 15, 2021, 3:47pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/12 "2021-07-15T15:47:22Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 17, 2021, 10:21pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/13 "2021-07-17T22:21:39Z")

</div>

Thanks! So I now have the following:

```julia
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?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 17, 2021, 11:17pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/14 "2021-07-17T23:17:20Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 18, 2021, 10:07pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/15 "2021-07-18T22:07:40Z")

</div>

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

```julia
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}

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 18, 2021, 10:13pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/16 "2021-07-18T22:13:45Z")

</div>

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

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 18, 2021, 10:56pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/17 "2021-07-18T22:56:47Z")

</div>

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

```julia
@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:

```julia
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](https://discourse.julialang.org/t/missing-from-the-variable-map-error/63172). Not sure if mine is exactly the same issue. I wasn’t able to figure it out…

---

<div class="post-metadata">

### Author: ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)
#### Post date: [July 19, 2021, 3:20am UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/18 "2021-07-19T03:20:05Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 19, 2021, 10:54pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/19 "2021-07-19T22:54:23Z")

</div>

Thanks! That was indeed the issue!

---

<div class="post-metadata">

### Author: ![Michael\_Barmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_barmann/32/20389_2.png) [@Michael\_Barmann](https://discourse.julialang.org/u/Michael_Barmann)
#### Post date: [July 19, 2021, 10:59pm UTC](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604/20 "2021-07-19T22:59:59Z")

</div>

Alright, here is my completed little toy problem!

```julia
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.

[Next page](https://discourse.julialang.org/t/stackoverflow-error-in-sir-model-using-modelingtookit/64604.md?page=2)