有限元方法教學探析

蘇海燕

（新疆大學 數學與系統科學學院，新疆烏魯木齊 830046）

在物理學、化學、生物學和數學等領域中，許多現象和規律通常可以用偏微分方程來揭示。然而，在一般情況下，絕大多數的偏微分方程初邊值問題無法給出解析解，采用有限元方法進行數值求解該類方程是當前主流的研究手段之一[1]。有限元課程的教學內容主要包括理論分析和數值模擬兩方面。在數值模擬方面，雖然現在已經有很多有限元數值模擬軟件 （如FreeFEM++, Fenics, COMSOL 等）僅僅利用偏微分方程所對應的變分形式，就可以通過幾行代碼非常便捷地得到數值解，但是這不利于初學者領悟有限元方法的核心和本質。

本文主要結合筆者的學習以及多年的有限元方法教學經驗，旨在通過引入方法，詳細地介紹有限元方法的主要思想，引導學生循序漸進地學習有限元方法的理論基礎和數值編程。

1 有限元方法

下面我們將詳細介紹有限元方法求解偏微分方程的主要步驟及其易于拓展到高維數和復雜問題的編程思想。

1.1 確定問題

本文考慮如下的橢圓形方程[2]：

1.2 變分格式

眾所周知，變分原理是有限元方法的基礎[3]，因此我們將首先推導(1)式的變分格式。 令

對（1）式左右兩端同時乘檢驗函數v，然后采用分部積分，則可得如下變分問題：

1.3 網格剖分

針對二維問題，網格剖分可分為三角形網格剖分、四邊形網格剖分、多邊形網格剖分等。采用不同的網格剖分方式，將直接影響下一節中基函數的建立。本文將采用最常見的三角形結構網格進行討論[4]。

1.4 基函數

由于方程(4)是無限維空間中的問題，有限元方法要求用有限維空間逼近無限維空間，在有限維空間中對原問題進行逼近[5]。 接下來介紹有限元空間的構造，假設

常見的多項式基函數有一次多項式基函數、 二次多項式基函數等。 在本文中，我們采用P1有限元，即基函數Φ 在每一片上是線性多項式Φ=ax+by+c，其中a，b，c 是常數。

1.5 離散格式

接下來，在上述建立的離散有限元空間中對連續變分格式（4）進行有限元離散[3]。 對，有

則可得離散格式，如下：

vh的任意性，取vh=Φj，j=1，2…，M，則有

綜上所述，問題（8）可以轉化為以下形式：

1.6 數值積分

由上一節的推導，可以發現(3)式是計算矩陣A，B，C，F 的基礎。 因此，如何快速、高精度地計算(3)式顯得尤為的重要。 接下來，這里簡單回顧一維高斯積分、二維正方形高斯積分和二維三角形高斯積分[6]。

若函數f（x）∈C2n[-1，1]，則使用高斯-勒讓德求積公式：

其中xk為高斯點，即是n 次勒讓德多項式Pn（x）的零點，An是求積系數：

若函數（x，y）∈C2n，2n[-1，1]2，則采用如下積分公式：

其中xi，xj分別是x，y 方向的高斯點，Ai，Aj為相應的求積系數。

若函數定義在單位三角形上時，需將積分進行坐標變換到正方形區域后利用(16)式進行計算：

其中

xi，xj分別是x，y 方向的高斯點，Ai，Aj為相應的求積系數。

1.7 組裝矩陣

在計算機編程實現中，數值計算矩陣A，B，C，F 通常有兩種思路： 第一種是對矩陣中的每個元素逐一進行計算，然后組裝為最后的總矩陣。第二種是對網格剖分中的每個單元進行矩陣的計算，最后進行組裝[7]。 此外，根據基函數的定義，我們可以發現在計算

的過程中，當Φi，Φj的非零函數區域不相交時有

因此，局部矩陣是稀疏的，為了節約計算儲存空間，我們在實際計算中利用小矩陣或者稀疏矩陣進行存儲，并記錄局部矩陣在總剛矩陣中的位置，為準確的總剛度矩陣合并提供保障。 另外矩陣AC 是對稱矩陣，所以可以通過矩陣的對稱性初步驗證矩陣組裝的正確性。 經過上述分塊組裝后，進行最后的合并：

1.8 邊界處理

在本文中，我們考慮第一類邊界條件，因此僅需對上一節中組裝好的總剛矩陣L 進行邊界處理。 需要將邊界節點i 所對應的行替換為0，再令對角元Iii=1，如：假設i 為邊界點，則進行如下處理：

此外，若邊界條件是第二類邊界條件時，在推導變分格式的過程中會有新的矩陣出現，然后將其新產生的矩陣組合在L，F 中即可，而不需要進行類似(22)式的處理過程。

1.9 求解Lx=b

經過之前的一系列處理，我們即將得到偏微分方程的數值解。接下來，對Lx=b 進行數值求解。針對矩陣的不同性質，采用不同的求解器。 本文旨在讓學生了解，采用有限元方法求解偏微分方程的具體步驟。 因此，直接調用程序庫中的一些求解器進行求解。常見的求解辦法有：Gauss 消去法，LU 分解法，Cholesky 分解法，迭代法，超松弛迭代法，共軛梯度法等[9]。

1.10 評估數值結果

雖然有了數值解，但我們對數值結果的優劣還未知。因此我們需要構造精確解算例，通過計算數值解的數值誤差來評估數值解的精確度。通常情況下，我們會計算如下三種誤差范數，以及結合網格參數計算收斂階，以此判斷數值解的準確性。

以及收斂階計算公式，

其中Ei，Ei+1 是對應于網格參數為hi，hi+1 時的誤差。 由文獻[10]可知當采用P1 元時，||u-uh||2L2≤C h2，|u-uh|H21≤C h 和||u-uh||H21C h，其中C 是給定的常數。

接下來，我們構造一個真解uexact=x（x-1）y（y-1），則f=x2y2+x2y+x2+xy2-3xy-x+y2-y。 我們對方程(1)進行有限元計算，結果如表1所示。

表1 有限元數值誤差結果及收斂階

由文獻[10]給出的誤差分析可知，我們的計算結果(見表1)符合理論分析。

2 結語

本文結合筆者的學習和教學經歷，對有限元方法的學習和教學過程進行了總結歸納。 以一個確定的偏微分方程，對如何設計編程實現有限元求解的每一個步驟進行了循序漸進的介紹。 雖然在編程上較為耗時且困難，但這個最經典的算法是大多數算法的基礎，需要初學者扎實掌握。