 # PDE Errors

I am trying to reproduce this solution to a PDE (Mathematica code located beneath tweet).

https://twitter.com/matthen2/status/1378820454589292550?s=21

This is my attempt thus far. Note that with my domains I actually have a rectangle (it seems a bit far off, but any help on how to make this an arbitrary shape such as an ellipse would also be fantastic).

``````using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

m = 5.0
@parameters t x y
@variables u(..)
Dt  = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = 1im * m * Dt(u(t,x,y)) ~ -Dxx(u(t, x, y)) - Dyy(u(t,x,y))

#bcs = [u(0, x, y) ~ exp(3im * (x+y) - 2 * (x^2+y^2)), Dirichlet0BC(Float64)]
bcs = vcat(u(0, x, y) ~ (3*cos(x+y) + 3im*sin(x+y)) * exp(-2.0 * (x^2 + y^2)), [Dirichlet0BC(Float64)])

domains = [t ∈ IntervalDomain(0.0, 15.0),
x ∈ IntervalDomain(0.0,  5.0),
y ∈ IntervalDomain(0.0,  3.0)]

pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])

dx = 0.01
dy = 0.01
order = 2
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)

prob = discretize(pdesys, discretization)
sol = solve(prob,Tsit5())
``````

The commented out boundary conditions are more along the lines of what I’d like to go for, but my mental model of the handling of complex numbers does not line up with that of Symbolics.jl (help on this would be great too!). The problem still seems to be with the definition of the boundary conditions (though the actual line with the error is the second from the bottom):

``````MethodError: no method matching operation(::Int64)
Closest candidates are:
operation(::Term) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:343
operation(::SymbolicUtils.Add) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:605
operation(::SymbolicUtils.Mul) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:748
...

Stacktrace:
 symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:80
 discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:152
 top-level scope
@ In:24
 eval
@ ./boot.jl:360 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094
``````

Notably, switching the order of the boundary conditions gets the error below, indicating that there’s probably issues with both boundary condition definitions.

``````type RobinBC has no field lhs
``````

I’m mostly basing this off of this example in the tests of DiffEqOperators.jl, which I was able to get working on my machine:

``````@parameters t x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
t_min= 0.
t_max = 2.0
x_min = 0.
x_max = 2.
y_min = 0.
y_max = 2.

# 3D PDE
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))

analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)
# Initial and boundary conditions
bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

# Space and time domains
domains = [t ∈ IntervalDomain(t_min,t_max),
x ∈ IntervalDomain(x_min,x_max),
y ∈ IntervalDomain(y_min,y_max)]
pdesys = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

# Method of lines discretization
dx = 0.1; dy = 0.2
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
prob = ModelingToolkit.discretize(pdesys,discretization)
sol = solve(prob,Tsit5())
``````

Any help would be greatly appreciated

1 Like

`Dirichlet0BC(Float64)` doesn’t make sense. What’s the condition supposed to be?

1 Like

Apologies, I am still learning both PDEs and the various tooling associated with them.

That was my attempt to translate this from the Mathematica code

``````DirichletCondition[u[t, xx, yy] == 0, True]
``````

My understanding was that the purpose of this is to make it zero at the spatial boundaries. Since I am using a rectangle, those could be put in directly, but eventually I’ll want to try different shapes so I was attempting to keep it general.

So then that would be a bunch of:

``````       u(t,x_min,y) ~ 0.0,
u(t,x_max,y) ~ 0.0,
u(t,x,y_min) ~ 0.0,
u(t,x,y_max) ~ 0.0
``````

as for making it general, right now you can’t. It’s one of the main things we are working on and you can see the suggested syntax in PDESystem boundary specification · Issue #526 · SciML/ModelingToolkit.jl · GitHub

1 Like

Awesome! That makes sense for the rectangular case. It seems generalizing to other shapes is making some substantial progress too. I’ve implemented this fix

``````m = 5.0
x_min, x_max = 0.0, 5.0
y_min, y_max = 0.0, 3.0

@parameters t x y
@variables u(..)
Dt  = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = 1im * m * Dt(u(t,x,y)) ~ -Dxx(u(t, x, y)) - Dyy(u(t,x,y))

bcs = vcat([u(t,x_min,y) ~ 0.0,
u(t,x_max,y) ~ 0.0,
u(t,x,y_min) ~ 0.0,
u(t,x,y_max) ~ 0.0],
u(0.0, x, y) ~ (3*cos(x+y) + 3im*sin(x+y)) * exp(-2.0 * (x^2 + y^2)))

domains = [t ∈ IntervalDomain(0.0, 15.0),
x ∈ IntervalDomain(x_min,  x_max),
y ∈ IntervalDomain(y_min,  y_max)]

pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])

dx = 0.01
dy = 0.01
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)

prob = discretize(pdesys, discretization)
sol = solve(prob,Tsit5())
``````

but the bug still persists

``````
MethodError: no method matching operation(::Int64)
Closest candidates are:
operation(::Term) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:343
operation(::SymbolicUtils.Add) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:605
operation(::SymbolicUtils.Mul) at /home/gmf/.julia/packages/SymbolicUtils/9ScfD/src/types.jl:748
...

Stacktrace:
 symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:80
 discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:152
 top-level scope
@ In:29
``````

The other bug appears to be gone, but it seems there’s something wrong with the t = 0 boundary condition. Thanks again!

Complex numbers are not supported in MOLFiniteDifference at this time. This is the underlying cause of your problem. Feel free to open an issue.

1 Like

Ah, that makes sense. Thanks for your help!