Hi everyone,
I’m trying out the new capability to deal with integrals in MethodOfLines.jl (in the new main, v0.7.6). I’m trying to simply expand an ODE example into a PDE. For those of you interested, the ODE is a simple SIR compartmental model, and the PDE is one where individuals are structured by age. It’s very close to working, except the sum of the states over time should remain constant, yet they decline over time (see the last line). I thought that setting the derivative wrt a to 0 at the upper boundary would prevent mass from leaving the system. Can anyone spot where I’m going wrong?
using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using Plots
β = 0.05
c = 10.0
γ = 0.25
S₀ = 990.0
I₀ = 10.0
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 26.0
amin = 0.0
amax = 52.0*75.0
@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
N(t) = Ia(S(a,t)) + Ia(I(a,t)) + Ia(R(a,t))
eqs = [Dt(S(a,t)) + Da(I(a,t)) ~ - β * c * S(a,t) * Ia(I(a,t))/N(t),
Dt(I(a,t)) + Da(I(a,t)) ~ β * c * S(a,t) * Ia(I(a,t))/N(t) - γ*I(a,t),
Dt(R(a,t)) + Da(R(a,t)) ~ γ*I(a,t)]
bcs = [
S(a,0.0) ~ S₀*da/(amax-amin),
S(0.0,t) ~ 0.0,
I(a,0.0) ~ I₀*da/(amax-amin),
I(0.0,t) ~ 0.0,
R(a,0.0) ~ 0.0,
R(0.0,t) ~ 0.0,
Da(S(amax,t)) ~ 0.0,
Da(I(amax,t)) ~ 0.0,
Da(R(amax,t)) ~ 0.0
]
domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]
@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t), R(a, t)])
discretization = MOLFiniteDifference([a=>da],t)
prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, Tsit5(), saveat = dt)
Smat = sol.u[S(a, t)]
Imat = sol.u[I(a, t)]
Rmat = sol.u[R(a, t)]
sum(Smat .+ Imat .+ Rmat, dims=1)