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:
[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.
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