(Experimental) Modeling Discrete Systems
In this example, we will use the new DiscreteSystem
API to create an SIR model.
using ModelingToolkit
using ModelingToolkit: t_nounits as t
using OrdinaryDiffEq: solve, FunctionMap
@inline function rate_to_proportion(r, t)
1 - exp(-r * t)
end
@parameters c δt β γ
@constants h = 1
@variables S(t) I(t) R(t)
k = ShiftIndex(t)
infection = rate_to_proportion(
β * c * I(k - 1) / (S(k - 1) * h + I(k - 1) + R(k - 1)), δt * h) * S(k - 1)
recovery = rate_to_proportion(γ * h, δt) * I(k - 1)
# Equations
eqs = [S(k) ~ S(k - 1) - infection * h,
I(k) ~ I(k - 1) + infection - recovery,
R(k) ~ R(k - 1) + recovery]
@mtkbuild sys = DiscreteSystem(eqs, t)
u0 = [S(k - 1) => 990.0, I(k - 1) => 10.0, R(k - 1) => 0.0]
p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1]
tspan = (0.0, 100.0)
prob = DiscreteProblem(sys, u0, tspan, p)
sol = solve(prob, FunctionMap())
retcode: Success
Interpolation: left-endpoint piecewise constant
t: 101-element Vector{Float64}:
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
⋮
92.0
93.0
94.0
95.0
96.0
97.0
98.0
99.0
100.0
u: 101-element Vector{Vector{Float64}}:
[0.24690087971667385, 989.5051237293776, 10.24797539090573]
[0.49992429364961877, 988.9982323978576, 10.501843308492768]
[0.7592157288009717, 988.479053494723, 10.761730776476044]
[1.0249238083995655, 987.9473092971556, 11.027766894444877]
[1.2972003431544508, 987.4027168203793, 11.300082836466316]
[1.5762003824739328, 986.8449877701768, 11.578811847349312]
[1.8620822655923719, 986.2738284979789, 11.864089236428692]
[2.1550076725435083, 985.6889399587333, 12.15605236872322]
[2.45514167491548, 985.090017671768, 12.454840653316563]
[2.7626527863189656, 984.4767516848771, 12.76059552880404]
⋮
[78.82720588986679, 843.9317615809414, 77.24103252919139]
[80.73429377803495, 840.6787392619304, 78.58696696003413]
[82.67461290570472, 837.3819011012713, 79.94348599302349]
[84.64842460763423, 834.0414203767291, 81.31015501563621]
[86.6559794879602, 830.6575018520915, 82.6865186599478]
[88.69751690774523, 827.2303822969169, 84.07210079533738]
[90.77326447234499, 823.7603309652128, 85.46640456244165]
[92.88343751961379, 820.2476500305196, 86.86891244986607]
[95.02823861000405, 816.692674974925, 88.27908641507038]
All shifts must be non-positive, i.e., discrete-time variables may only be indexed at index k, k-1, k-2, ...
. If default values are provided, they are treated as the value of the variable at the previous timestep. For example, consider the following system to generate the Fibonacci series:
@variables x(t) = 1.0
@mtkbuild sys = DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t)
\[ \begin{align} Shift(t, 1)\left( \mathtt{x_{t - 1}}\left( t \right) \right) &= x\left( t \right) \\ Shift(t, 1)\left( x\left( t \right) \right) &= \mathtt{x_{t - 1}}\left( t \right) + x\left( t \right) \end{align} \]
The "default value" here should be interpreted as the value of x
at all past timesteps. For example, here x(k-1)
and x(k-2)
will be 1.0
, and the initial value of x(k)
will thus be 2.0
. During problem construction, the past value of a variable should be provided. For example, providing [x => 1.0]
while constructing this problem will error. Provide [x(k-1) => 1.0]
instead. Note that values provided during problem construction do not apply to the entire history. Hence, if [x(k-1) => 2.0]
is provided, the value of x(k-2)
will still be 1.0
.