Please try to make a minimum working example by stripping out as much as possible while keeping the same error, that will make it much easier to help
I stripped away some unnecessary stuff. The method works if I have only t and a as variables. But in this version there is another continous variable b.
This is the code now:
using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
tmin = 0.0
tmax = 50.0
amin = 0.0
amax = 100.0
bmin = 0.0
bmax = 1.0
#Model parameters
β = 3.0 #transmission rate
γ = 1.0 # recovery rate
alpha = 0.05 #waning rate
N = 100.0 #population size
k = 10.0
lambda = 90.0
mu = 0.01
#Heaviside function
theta(a) = (a > 0) ? 1 : 0
@register_symbolic theta(a)
# Parameters, variables, and derivatives
@parameters t a b
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Db = Differential(b)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
Ib = Integral(b in DomainSets.ClosedInterval(bmin,bmax))
# PDE and boundary conditions
eqs = [
Dt(S(t,a,b)) + Da(S(t,a,b)) + Db(S(t,a,b)) ~ - β * S(t,a,b) * Ia(Ib(I(t,a,b)))/N - mu*S(t,a,b),
Dt(I(t,a,b)) + Da(I(t,a,b)) + Db(I(t,a,b)) ~ β * S(t,a,b) * Ia(Ib(I(t,a,b)))/N - γ*I(t,a,b) - mu*I(t,a,b),
Dt(R(t,a,b)) + Da(R(t,a,b)) + Db(R(t,a,b)) ~ γ * I(t,a,b) - mu*R(t,a,b)]
S₀ = N * 0.99 #initial number susceptible
I₀ = N * 0.01 #initial number infected
R₀ = 0.0 #initial number recovered
min_age = 0.0 #minimum age
max_age = 100.0 #maximum age
min_b = 0.0 #minimum frailty
max_b = 1.0 #maximum frailty
bcs = [
S(0.0,a,b) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b),
I(0.0,a,b) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b),
R(0.0,a,b) ~ 0.0,
S(t,0.0,b) ~ S₀/(max_age - min_age),
I(t,0.0,b) ~ I₀/(max_age - min_age),
R(t,0.0,b) ~ 0.0,
S(t,a,0.0) ~ N/(max_age - min_age),
I(t,a,0.0) ~ 0.0,
R(t,a,0.0) ~ 0.0
]
# Space, age, and frailty domains
domains = [t ∈ Interval(tmin,tmax),a ∈ Interval(amin,amax), b ∈ Interval(bmin,bmax)]
# PDE system
@named SIRmodel = PDESystem(eqs, bcs, domains, [t,a,b], [S(t,a,b), I(t,a,b), R(t,a,b)])
# Method of lines discretization
da = 10.0
db = 0.1
discretization = MOLFiniteDifference([a=>da,b=>db],t)
# Convert the PDE problem into an ODE problem
prob = discretize(SIRmodel, discretization)
# Solve ODE problem
dt = tmax/100.0
sol_pde = solve(prob, Tsit5(), saveat = dt)
When I try to run you code on Julia v1.9.3 and MethodOfLines v0.9.6, I first get a different error:
prob = discretize(SIRmodel, discretization)
ERROR: AssertionError: Integral Integral(a, 0.0 .. 100.0)(Integral(b, 0.0 .. 1.0)(I(t, a, b))) must be purely of a variable, got Integral(b, 0.0 .. 1.0)(I(t, a, b)). Try wrapping the integral argument with an auxiliary variable.
I think the issue is explained here:
https://docs.sciml.ai/MethodOfLines/stable/tutorials/PIDE/#integral
Modifying the code as below solves that error (but not the original one you mentioned). Note the two new auxiliary fields which represent the integrals.
Fix for integrals
using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
tmin = 0.0
tmax = 50.0
amin = 0.0
amax = 100.0
bmin = 0.0
bmax = 1.0
#Model parameters
β = 3.0 #transmission rate
γ = 1.0 # recovery rate
alpha = 0.05 #waning rate
N = 100.0 #population size
k = 10.0
lambda = 90.0
mu = 0.01
#Heaviside function
theta(a) = (a > 0) ? 1 : 0
@register_symbolic theta(a)
# Parameters, variables, and derivatives
@parameters t a b
@variables S(..) I(..) R(..) ItotalA(..) Itotal(..)
Dt = Differential(t)
Da = Differential(a)
Db = Differential(b)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
Ib = Integral(b in DomainSets.ClosedInterval(bmin,bmax))
# PDE and boundary conditions
eqs = [
ItotalA(t,a,b) ~ Ib(I(t,a,b)),
Itotal(t,a,b) ~ Ia(ItotalA(t,a,b)),
Dt(S(t,a,b)) + Da(S(t,a,b)) + Db(S(t,a,b)) ~ - β * S(t,a,b) * Itotal(t,a,b)/N - mu*S(t,a,b),
Dt(I(t,a,b)) + Da(I(t,a,b)) + Db(I(t,a,b)) ~ β * S(t,a,b) * Itotal(t,a,b)/N - γ*I(t,a,b) - mu*I(t,a,b),
Dt(R(t,a,b)) + Da(R(t,a,b)) + Db(R(t,a,b)) ~ γ * I(t,a,b) - mu*R(t,a,b)]
S₀ = N * 0.99 #initial number susceptible
I₀ = N * 0.01 #initial number infected
R₀ = 0.0 #initial number recovered
min_age = 0.0 #minimum age
max_age = 100.0 #maximum age
min_b = 0.0 #minimum frailty
max_b = 1.0 #maximum frailty
bcs = [
S(0.0,a,b) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b),
I(0.0,a,b) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b),
R(0.0,a,b) ~ 0.0,
S(t,0.0,b) ~ S₀/(max_age - min_age),
I(t,0.0,b) ~ I₀/(max_age - min_age),
R(t,0.0,b) ~ 0.0,
S(t,a,0.0) ~ N/(max_age - min_age),
I(t,a,0.0) ~ 0.0,
R(t,a,0.0) ~ 0.0
]
# Space, age, and frailty domains
domains = [t ∈ Interval(tmin,tmax),a ∈ Interval(amin,amax), b ∈ Interval(bmin,bmax)]
# PDE system
@named SIRmodel = PDESystem(eqs, bcs, domains, [t,a,b], [S(t,a,b), I(t,a,b), R(t,a,b), ItotalA(t,a,b), Itotal(t,a,b)])
# Method of lines discretization
da = 10.0
db = 0.1
discretization = MOLFiniteDifference([a=>da,b=>db],t)
# Convert the PDE problem into an ODE problem
prob = discretize(SIRmodel, discretization)
# Solve ODE problem
dt = tmax/100.0
sol_pde = solve(prob, Tsit5(), saveat = dt)
From looking at the stacktrace, it seems like the error appears when evaluating the ODE function. Running just the function gives the same error:
prob.f(copy(prob.u0), prob.u0, Float64[], 0.0)
ERROR: MethodError: Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64} to an object of type Float64
I think the problem is indeed with the double integrals. Here is a slightly shortened MWE to reproduce the error. The symbolically discretized ODEs still contain the symbolic Integral expressions.
Reduced MWE
using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
tmin = 0.0
tmax = 50.0
amin = 0.0
amax = 100.0
bmin = 0.0
bmax = 1.0
N = 100.0
S₀ = N
@parameters t a b
@variables S(..) StotalA(..) Stotal(..)
Dt = Differential(t)
Da = Differential(a)
Db = Differential(b)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
Ib = Integral(b in DomainSets.ClosedInterval(bmin,bmax))
# Want to use integral of S over whole a-b-domain, using two auxiliary variables.
eqs = [
StotalA(t,a,b) ~ Ib(S(t,a,b)),
Stotal(t,a,b) ~ Ia(StotalA(t,a,b)),
Dt(S(t,a,b)) + Da(S(t,a,b)) + Db(S(t,a,b)) ~ - S(t,a,b) * Stotal(t,a,b) / N - S(t,a,b),
]
bcs = [
S(0.0,a,b) ~ S₀,
S(t,0.0,b) ~ S₀,
S(t,a,0.0) ~ N,
]
domains = [t ∈ Interval(tmin,tmax), a ∈ Interval(amin,amax), b ∈ Interval(bmin,bmax)]
@named SIRmodel = PDESystem(eqs, bcs, domains, [t,a,b], [S(t,a,b), StotalA(t,a,b), Stotal(t,a,b)])
discretization = MOLFiniteDifference([a=>amax/10,b=>bmax/10],t)
prob = discretize(SIRmodel, discretization)
# Produces the error
# prob.f(copy(prob.u0), prob.u0, Float64[], 0.0)
probS = symbolic_discretize(SIRmodel, discretization)
# Integral along b axis is properly discretized
@show probS[1].eqs[1]
# (probS[1]).eqs[1] = (StotalA(t))[1, 1] - 0.05(S(t))[1, 1] - 0.1(S(t))[1, 2] - 0.1(S(t))[1, 3] - 0.1(S(t))[1, 4] - 0.1(S(t))[1, 5] - 0.1(S(t))[1, 6] - 0.1(S(t))[1, 7] - 0.1(S(t))[1, 8] - 0.1(S(t))[1, 9] - 0.1(S(t))[1, 10] - 0.05(S(t))[1, 11] ~ 0
# Integral along a axis is not discretized
@show probS[1].eqs[122]
# (probS[1]).eqs[122] = (Stotal(t))[1, 1] - Integral(a, 0.0 .. 100.0)((StotalA(t))[1, 1]) ~ 0
Should double integrals be working at this point? Removing the second integration seems to discretize and solve just fine.
Thank you very much, Sevi! I will look into these issues and your suggestions in the coming days.