【我已經發了同名西瓜視頻，本文僅是視頻中Markdown文件内容，不包括視頻中講解】

如果你希望和我一樣Win10上安裝FEniCS，建議采用Win10 WSL2 Ubuntu。

FEniCS是一個有限元法求解偏微分方程的開源計算平台。

• 泊松方程
• 将泊松方程轉化成（弱）變分形式
• 求解微分方程邊值問題（1維泊松方程）
• 求解2維泊松方程
• ParaView中查看解

泊松方程

第一個例子自然是求解著名的泊松方程。

将泊松方程轉化成（弱）變分形式

首先：對泊松方程兩邊都乘上測試函數，然後對全域​積分。

這裡需要說明一下試探函數和測試函數的概念。

測試函數，定義為：

試探函數，定義為：

然後，對上面這個變分方程進行變換

左邊是雙線性形式​，右邊是線性形式​:

于是得到标準的變分問題：尋求​，滿足

這就是和前面的泊松問題等價的變分問題。

求解微分方程邊值問題（1維泊松方程）

有了前面的準備，現在可以寫代碼求解了。

第一步：為求解域創建網格，并定義函數空間

from dolfin import * # 單元個數 nel = 20 # 左右端點 xmin = 0 xmax = 1 # 試探/測試函數的多項式階數 p = 2 # 創建網格 mesh = IntervalMesh(nel,xmin,xmax) # 使用連續Galerkin定義函數空間 # 每個單元上的p階(拉格朗日)函數 V = FunctionSpace(mesh,"CG",p)

第二步：定義已提供的表達式

u0 = Constant(0) f = Constant(-1) g = Constant(1)

第三步：創建和應用Dirichlet邊界條件

# 定義Dirichlet邊界(x = 0) def dirichlet_boundary(x, on_boundary): return on_boundary and abs(x[0]) < DOLFIN_EPS # 在點x=0施加一個Dirichlet條件 bc = DirichletBC(V,u0,dirichlet_boundary)

第四步：定義變分問題

u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u),grad(v))*dx L = f*v*dx g*v*ds

第五步：求解并繪圖

# 求解 u = Function(V) solve(a == L, u, bc) # 繪圖 plot(u)

求解2維泊松方程

• ​ (一個單位矩陣)
• ​(Dirichlet邊界)
• ​ (Neumann邊界)
• ​ (法向導數)
• ​ (源項)

類似的，也分五步求解

from dolfin import * # 第一步：為求解域創建網格，并定義函數空間 mesh = UnitSquareMesh(32, 32) V = FunctionSpace(mesh, "Lagrange", 1) # 第二步：定義已提供的表達式 u0 = Constant(0.0) f = Expression("10*exp(-(pow(x[0] - 0.5, 2) pow(x[1] - 0.5, 2)) / 0.02)", degree=2) g = Expression("sin(5*x[0])", degree=2) # 第三步：創建和應用`Dirichlet`邊界條件 def boundary(x): return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS bc = DirichletBC(V, u0, boundary) # 第四步：定義變分問題 u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u), grad(v))*dx L = f*v*dx g*v*ds # 第五步：求解并繪圖 u = Function(V) solve(a == L, u, bc) plot(u)

ParaView中查看解

# 将解保存為VTK格式 file = File("poisson.pvd") file << u

【結束】

,

更多精彩资讯请关注tft每日頭條，我们将持续为您更新最新资讯!

查看全部