華中師范大學(xué) hahakity

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

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

學(xué)習(xí)目標(biāo)Galerkin 有限元法

分片基底函數(shù)展開

二階微分算子的降階

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

預(yù)備知識加權(quán)殘差法 (Weighted Residual Method)

Galerkin 有限元法

假設(shè)要求解的微分方程有統(tǒng)一的格式，

其中

是微分算子，

是要求解的場，

是源。

將

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

求和符號)，

注意此時

為基底函數(shù)，人工設(shè)定，為已知量，

為待定系數(shù)。

選取基底函數(shù)

作為伽遼金法中的檢驗函數(shù)，則對第 i 個檢驗函數(shù)，加權(quán)殘差為，

其中

是全域，

,

可以將

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

表示列向量，

是

構(gòu)成的矩陣，

是

構(gòu)成的列向量，

是

構(gòu)成的列向量。

Galerkin 有限元法令殘差

為0，通過解

，得到待定系數(shù)

。

分片基底函數(shù)展開

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

近似為有限個單元上的分片函數(shù)。如果單元很小，每個單元內(nèi)

近似線性變化。

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

。

對第 e 個單元區(qū)域

，線性插值函數(shù)為

其中

，

當(dāng)

時，

,

當(dāng)

時，

,

根據(jù)

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

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

在第 j 個單元有重疊，跟

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

根據(jù)

可知第 j 行只有

,

,

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

二階微分算子降階

算

時，很快遇到第一個問題，線性插值函數(shù)一階微分在每個單元內(nèi)為常數(shù)，

二階或更高階導(dǎo)數(shù)等于0。如果微分算子里有二階導(dǎo)數(shù)，可以用分部積分法降階，

二階微分算子降階之后變?yōu)?#xff0c;

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

或它的導(dǎo)數(shù)已知，因此可以將這一項單獨列出來，令

在新的定義下，方程變?yōu)?#xff0c;

多出來的那項

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

構(gòu)成的列向量。

此時矩陣

第 j 行的非零矩陣元為,

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

注意此時

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

簡單的泊松方程求解

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

選擇基底函數(shù)，對每個子域插值

用 Galerkin 法改寫微分方程，函數(shù)內(nèi)積計算

的矩陣元

解

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

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

，而電場本身又可寫成標(biāo)量勢

的負(fù)梯度

, 其中

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

，二維情況下方程為，

，

一維情況下

其中

是真空電極化常數(shù)。

選擇求解區(qū)域

，邊界條件

。

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

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

, 我們先假裝不知道。

計算

的矩陣元，

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

與

。在 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)

輸出結(jié)果：

Jackson 電動力學(xué)書本中介紹了一種簡化方法，假設(shè)

在每個單元里是常數(shù)，只對

積分，則積分結(jié)果為，

先不管

, 矩陣方程的解為，

。寫一段簡單的 python 代碼，構(gòu)造矩陣

和列向量

，使用 scipy.sparse.linalg.inv 函數(shù)求

的逆，進而得到結(jié)果。

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 稀疏矩陣庫，構(gòu)造 3 對角稀疏矩陣 K

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

return K

def bvec(self):

'''假設(shè) rho(x) 在每個單元內(nèi)值為常數(shù)，僅對 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 個節(jié)點，計算結(jié)果比較精確，101 或 11 個節(jié)點，計算結(jié)果誤差較大。

fem = FEM(1001)

fem.compare()

二維或高維時的網(wǎng)格生成

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

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

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

pip install gmsh

pip install pygmsh

使用 pygmsh，可以很快生成一個二維的圓形區(qū)域網(wǎng)格，

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()

# 將網(wǎng)格文件保存為 vtk 格式

mesh.write("test.vtk")

安裝 paraview，打開 test.vtk, 則可以看到如下二維圓形區(qū)域的網(wǎng)格圖，

總結(jié)

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

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

參考文獻：

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。