華中師范大學 hahakity

第一部分介紹了加權殘差法(Weighted Residual Method)并引出有限元算法中經常出現的伽遼金(Galerkin)法。簡單說明了有限元與有限差分法的區別：有限差分是對微分算子做差分近似，有限元是對待求函數做基函數展開近似。這一節繼續介紹伽遼金有限元算法。

這里特意強調伽遼金有限元算法，是因為還有另一種基于泛函和變分的 Ritz 有限元算法。跳過 Ritz 變分法，直接講解伽遼金法入門比較簡單。

學習目標Galerkin 有限元法

分片基底函數展開

二階微分算子的降階

一維泊松方程：已知電荷密度求電勢

預備知識加權殘差法 (Weighted Residual Method)

Galerkin 有限元法

假設要求解的微分方程有統一的格式，

其中

是微分算子，

是要求解的場，

是源。

將

做函數基底展開 (使用愛因斯坦求和規則，省略

求和符號)，

注意此時

為基底函數，人工設定，為已知量，

為待定系數。

選取基底函數

作為伽遼金法中的檢驗函數，則對第 i 個檢驗函數，加權殘差為，

其中

是全域，

,

可以將

寫成矩陣乘積形式，此處用

表示列向量，

是

構成的矩陣，

是

構成的列向量，

是

構成的列向量。

Galerkin 有限元法令殘差

為0，通過解

，得到待定系數

。

分片基底函數展開

將全域離散化為 n 個子域，每個子域稱作一個單元(“element”)。連續的函數

近似為有限個單元上的分片函數。如果單元很小，每個單元內

近似線性變化。

先考慮簡單的一維問題，將區間分成 n 等份，每份長

。

對第 e 個單元區域

，線性插值函數為

其中

，

當

時，

,

當

時，

,

根據

的定義，得到第 j 個基底函數(注意它橫跨兩個單元)，

的函數圖像就是下圖中間那個大三角形，覆蓋第 j 和第 j+1 個單元。它跟

在第 j 個單元有重疊，跟

在第 j+1 個單元有重疊。此時內積都是局部進行，不重疊區域積分為0。

根據

可知第 j 行只有

,

,

三個非零矩陣元，K 矩陣是稀疏的 3 對角陣。

二階微分算子降階

算

時，很快遇到第一個問題，線性插值函數一階微分在每個單元內為常數，

二階或更高階導數等于0。如果微分算子里有二階導數，可以用分部積分法降階，

二階微分算子降階之后變為，

只與全域邊界條件有關，對于大部分不在邊界上的點，這項為0。剛好在邊界上的那些點，

或它的導數已知，因此可以將這一項單獨列出來，令

在新的定義下，方程變為，

多出來的那項

是分部積分項在邊界處的貢獻

構成的列向量。

此時矩陣

第 j 行的非零矩陣元為,

每個矩陣元的積分區間由前面的示意圖給出。

注意此時

或者等于 1/h，或者等于 -1/h，積分很好算。

簡單的泊松方程求解

我們已經有了足夠的知識，從頭構造一個簡單的一維 Galerkin 有限元算法。這里先演示一維有限元算法。步驟總結如下，將全域離散化為有限個子域(又稱單元)

選擇基底函數，對每個子域插值

用 Galerkin 法改寫微分方程，函數內積計算

的矩陣元

解

這里拿一個簡單的一維泊松方程演示伽遼金有限元算法。

物理問題：根據 Maxwell 方程，電場的散度正比于電荷密度

，而電場本身又可寫成標量勢

的負梯度

, 其中

。將第二個方程代入第一個得到給定電荷分布下靜電勢滿足的柏松方程，

，二維情況下方程為，

，

一維情況下

其中

是真空電極化常數。

選擇求解區域

，邊界條件

。

選擇一個簡單的電荷密度分布，

此時，泊松方程有簡單的解析解：

, 我們先假裝不知道。

計算

的矩陣元，

這個積分用手積挺容易出錯，寫一小段 python 代碼，將積分區域劃分為

與

。在 jupyter-notebook 里運行，

import sympy

xjm1, xj, xjp1 = sympy.symbols(['x_{j-1}', 'x_{j}', 'x_{j+1}'])

x, h, eps = sympy.symbols(['x', 'h', '\epsilon'])

rho = lambda x: sympy.sin(sympy.pi * x)

A = - sympy.integrate(rho(x) * (x - xjm1)/h, (x, xjm1, xj))

B = - sympy.integrate(rho(x) * (xjp1 - x)/h, (x, xj, xjp1))

sympy.simplify(A + B)

輸出結果：

Jackson 電動力學書本中介紹了一種簡化方法，假設

在每個單元里是常數，只對

積分，則積分結果為，

先不管

, 矩陣方程的解為，

。寫一段簡單的 python 代碼，構造矩陣

和列向量

，使用 scipy.sparse.linalg.inv 函數求

的逆，進而得到結果。

from scipy.sparse import dia_matrix

from scipy.sparse.linalg import inv

from numpy import pi

class FEM:

def __init__(self, nodes, xmin=0, xmax=1):

self.nodes = nodes

x = np.linspace(xmin, xmax, nodes)

self.x = x

self.h = x[1] - x[0]

def Kmatrix(self):

n = self.nodes

m = 1/self.h * np.ones(n)

data = [m, -2*m, m]

offsets = [-1, 0, 1]

# 使用 scipy.sparse 稀疏矩陣庫，構造 3 對角稀疏矩陣 K

K = dia_matrix((data, offsets), shape=(n, n)).tocsc()

return K

def bvec(self):

'''假設 rho(x) 在每個單元內值為常數，僅對 u_j(x) 做積分'''

return - np.sin(pi * self.x) * self.h

def solve(self):

K = self.Kmatrix()

b = self.bvec()

return inv(K) * b

def compare(self):

ground_truth = 1/(pi**2) * np.sin(pi * self.x)

fem_res = self.solve()

plt.plot(self.x, ground_truth, label="ground truth")

plt.plot(self.x, fem_res, label="finite element")

plt.title("number of nodes =%s"%self.nodes)

plt.xlabel("x")

plt.ylabel(r"$\phi(x)$")

plt.legend(loc='best')

plt.show()

如果用 1001 個節點，計算結果比較精確，101 或 11 個節點，計算結果誤差較大。

fem = FEM(1001)

fem.compare()

二維或高維時的網格生成

上面的一維問題有限元與有限差分算法區別很小，將方程右邊的

除到左邊，則 K 矩陣變成了二階導數的差分矩陣。但對于高維且有復雜邊界的偏微分問題，有限元算法就顯示出它的優勢。這一篇變得太長，下一篇嘗試構造通用的二維和三維有限元求解器，以及如何在二維圓形區域求解本文例子中的靜電勢。

為了簡化任務，使用 gmsh 軟件的 python 庫生成有限元分析需要的二維和三維網格。安裝方法如下，

pip install gmsh

pip install pygmsh

使用 pygmsh，可以很快生成一個二維的圓形區域網格，

import pygmsh

with pygmsh.geo.Geometry() as geom:

geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)

mesh = geom.generate_mesh()

# 將網格文件保存為 vtk 格式

mesh.write("test.vtk")

安裝 paraview，打開 test.vtk, 則可以看到如下二維圓形區域的網格圖，

總結

介紹了伽遼金有限元算法，一維分片函數展開，二階微分算子降階，數值求解了簡單的一維泊松方程，介紹了網格生成軟件 gmsh。將二維和三維有限元算法的通用求解器推遲到下一節。

警告：本文的求解器非通用求解器，僅作原理展示使用。

參考文獻：

【1】Tao Pang，An Introduction to computational physics (second edition)

總結

以上是生活随笔為你收集整理的python简单的计算方法_用python实现简单的有限元方法（二）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。