# PDE discretization error in MethodOfLines.jl

**URL:** <https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566>\
**Category:** Modelling & Simulations\
**Created:** [March 27, 2022, 6:53pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566 "2022-03-27T18:53:50Z")\
**Posts on this page:** 4\
**Page:** 2

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [September 4, 2023, 3:18pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/21 "2023-09-04T15:18:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [September 5, 2023, 12:57pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/22 "2023-09-05T12:57:49Z")

</div>

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)

```

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [September 13, 2023, 2:48pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/23 "2023-09-13T14:48:33Z")

</div>

When I try to run you code on Julia v1.9.3 and MethodOfLines v0.9.6, I first get a different error:

```julia
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](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**
>
> ```julia
> 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:

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

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [September 15, 2023, 11:27am UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/24 "2023-09-15T11:27:58Z")

</div>

Thank you very much, Sevi! I will look into these issues and your suggestions in the coming days.

[Previous page](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566.md?page=1)