# PDE Errors
**URL:** https://discourse.julialang.org/t/pde-errors/59707
**Category:** Modelling & Simulations
**Tags:** question, diffeq, pde, ode, symbolic
**Created:** [April 21, 2021, 1:41am UTC](https://discourse.julialang.org/t/pde-errors/59707 "2021-04-21T01:41:22Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)
#### Post date: [April 21, 2021, 1:41am UTC](https://discourse.julialang.org/t/pde-errors/59707/1 "2021-04-21T01:41:22Z")
I am trying to reproduce this solution to a PDE (Mathematica code located beneath tweet).
[https://twitter.com/matthen2/status/1378820454589292550?s=21](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).
```julia
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):
```julia
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:
[1] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:80
[2] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:152
[3] top-level scope
@ In[175]:24
[4] eval
@ ./boot.jl:360 [inlined]
[5] 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.
```julia
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:
```julia
@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
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 21, 2021, 10:46am UTC](https://discourse.julialang.org/t/pde-errors/59707/2 "2021-04-21T10:46:15Z")
> [@Ian\_Slagle](#):
>
> ```julia
> bcs = vcat(u(0, x, y) ~ (3*cos(x+y) + 3im*sin(x+y)) * exp(-2.0 * (x^2 + y^2)), [Dirichlet0BC(Float64)])
>
> ```
`Dirichlet0BC(Float64)` doesn’t make sense. What’s the condition supposed to be?
---
### Author:  [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)
#### Post date: [April 21, 2021, 1:34pm UTC](https://discourse.julialang.org/t/pde-errors/59707/3 "2021-04-21T13:34:30Z")
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
```julia
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 21, 2021, 7:19pm UTC](https://discourse.julialang.org/t/pde-errors/59707/4 "2021-04-21T19:19:30Z")
So then that would be a bunch of:
```julia
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](https://github.com/SciML/ModelingToolkit.jl/issues/526)
---
### Author:  [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)
#### Post date: [April 22, 2021, 1:08am UTC](https://discourse.julialang.org/t/pde-errors/59707/5 "2021-04-22T01:08:51Z")
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
```julia
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
```julia
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:
[1] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:80
[2] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/Fmo6r/src/MOLFiniteDifference/MOL_discretization.jl:152
[3] top-level scope
@ In[7]:29
```
The other bug appears to be gone, but it seems there’s something wrong with the t = 0 boundary condition. Thanks again!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 22, 2021, 6:02am UTC](https://discourse.julialang.org/t/pde-errors/59707/6 "2021-04-22T06:02:43Z")
Complex numbers are not supported in MOLFiniteDifference at this time. This is the underlying cause of your problem. Feel free to open an issue.
---
### Author:  [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)
#### Post date: [April 22, 2021, 4:36pm UTC](https://discourse.julialang.org/t/pde-errors/59707/7 "2021-04-22T16:36:34Z")
Ah, that makes sense. Thanks for your help!