我想尝试用Julia求解偏微分方程,有限差分法;
例子是在这里找到的: https://diffeqoperators.sciml.ai/stable/symbolic_tutorials/mol_heat/
我想把原文中的矩形域换成其他的形式,比如圆形域,因为我是初学者,搞来搞去也不知道怎么弄,从官方的文档里面也没有找到有用的信息?
域的定义和DomainSets.jl应该怎样设置?有Julia处理偏微分方程相关的教程或资料推荐吗?
我想读这个包的源码,查看Julia包源码的正确姿势应该是什么?谢谢诸位!
# Initial
domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0.0..1.0), Symbolics.VarDomainPairing(y, 0.0..1.0)]
domains = [(x, y) ∈ UnitDisk()]
@parameters x y
domains = [x ∈ UnitCircle(), y ∈ UnitCircle()]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, UnitCircle()), Symbolics.VarDomainPairing(y, UnitCircle())]
domains = [(x, y) ∈ (0 .. 1) × (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing((x, y), (0 .. 1) × (0 .. 1))]
domains = [x ∈ (0 .. 1) × (0 .. 1), y ∈ (0 .. 1) × (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, (0 .. 1) × (0 .. 1)), Symbolics.VarDomainPairing(y, (0 .. 1) × (0 .. 1))]
# it works!
domains = [x ∈ (0 .. 1), y ∈ (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0..1), Symbolics.VarDomainPairing(y, 0..1)]
using ModelingToolkit, DiffEqOperators, LinearAlgebra, DomainSets
using ModelingToolkit: Differential
using DifferentialEquations
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
dx = 0.1
dy = 0.1
bcs = [u(0, y) ~ 0.0,
u(1, y) ~ y,
u(x, 0) ~ 0.0,
u(x, 1) ~ x]
# Space and time domains
domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])
# Note that we pass in `nothing` for the time variable `t` here since we
# are creating a stationary problem without a dependence on time, only space.
discretization = MOLFiniteDifference([x => dx, y => dy], nothing, centered_order = 2)
prob = discretize(pdesys, discretization)
sol = solve(prob)
using Plots
xs, ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
u_sol = reshape(sol.u, (length(xs), length(ys)))
plot(xs, ys, u_sol, linetype = :contourf, title = "solution")