有限元數值解法在MATLAB中的實現及可視化

摘 要：偏微分方程的數值解法在數值分析中占有很重要的地位，很多科學技術問題的數值計算包括了偏微分方程的數值解問題。在學習初等函數時，總是先畫出它們的圖形，因為圖形能幫助我們了解函數的性質。而對于偏微分方程，畫出它們的圖形并不容易，尤其是沒有解析解的偏微分方程，畫圖就顯得更加不容易了。為了從偏微分方程的數學表達式中看出其所表達的圖形、函數值與自變量之間的關系，通過MATLAB編程，用有限元數值解法求解了偏微分方程，并將其結果可視化。

關鍵詞：偏微分方程；MATLAB；有限元法；可視化

中圖分類號：TP311.12 文獻標識碼：A

1 引言（Introduction）

偏微分方程的數值解法在數值分析中占有很重要的地位，很多科學技術問題的數值計算包括了偏微分方程的數值解問題。近三十多年來，它的理論和方法都有了很大的發展，而且在各個科學技術的領域中應用也愈來愈廣泛。例如，核武器的研制要有理論設計和核試驗。但核反應和核爆炸的過程是在高溫高壓的條件下進行的，而且巨大的能量在極短的時間內釋放出來，核裝置內部的細致反應過程及各個物理量的變化是根本不能用儀器測量出來的，核試驗只是提供綜合的數據。而描述核反應和爆炸物理過程的數學模型是一個很復雜的非線性偏微分方程組，也根本沒有辦法得到這個方程組理論上的精確解。所以發展核武器的國家都在計算機上對核反應過程進行數值模擬，這也稱為“數值核實驗”，它可以大大減少核試驗的次數，節約大量的經費，縮短研制的周期[1]。

在學習初等函數時，總是先畫出它們的圖形，因為圖形能幫助我們了解函數的性質。而對于偏微分方程，畫出它們的圖形并不容易，尤其是沒有解析解的偏微分方程，畫圖就顯得更加不容易了。所以本文主要研究如何用MATLAB數值求解偏微分方程，并將其數值解繪制成三維圖形的形式，從而可以從復雜的數學表達式中看出其所表達的圖像、函數值與自變量之間的關系[2]。

2 有限元法（Finite element method）

2.1 有限元法概述

有限元法的基本思想是將結構離散化，用有限個容易分析的單元來表示復雜的對象，單元之間通過有限個節點相互連接，然后根據變形協調條件綜合求解。由于單元的數目是有限的，節點的數目也是有限的，所以稱為有限元法。

一般來說，用差分法解偏微分方程，解得的結果就是方程的準確解函數在節點上的近似值。而用變分近似方法求解，是將近似解表示成有限維子空間中基函數的線性組合。在古典變分方法中，這樣的基函數一般采用冪函數和三角函數等初等函數，又要求在區域的邊界上滿足邊界條件，如果是二維或三維的不規則區域，這樣的基函數往往很難構造出來。所以，古典的變分方法雖然是得到近似的解析解（與差分方法不同），但是對一般的區域，卻往往難以實現。有限元方法，也是基于變分原理，由于選擇了特殊的基函數，使它能適用于較一般的區域。這種基函數是區域的剖分有關的，近似解u表示為基函數的線性組合，而線性組合中的系數，又是剖分節點上u或其導數的近似解。所以有限元方法既是基于變分原理，又具有差分方法的一些特點，并且適合于較復雜的區域和不同粗細的網格。正是由于這些特點，20世紀60年代以來，有限元方法的理論和應用得到迅速的發展，適用范圍也愈來愈廣泛。

2.2 有限元法的基本思想

有限元法是用有限個單元將連續體離散化，通過對有限個單元作分片插值求解各種力學、物理問題的一種數值方法。有限元法把連續體離散成有限個單元：桿系結構的單元是每一個桿件；連續體的單元是各種形狀（如三角形、四邊形、六面體等）的單元體。每個單元的場函數是只包含有限個待定節點參量的簡單場函數，這些單元場函數的集合就能近似代表整個連續體的場函數。根據能量方程或加權殘量方程可建立有限個待定參量的代數方程組，求解此離散方程組就得到有限元法的數值解。

有限元方法的基礎是變分原理和加權余量法，其基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內，選擇一些合適的節點作為求解函數的插值點，將微分方程中的變量改寫成由各變量或其導數的節點值與所選用的插值函數組成的線性表達式，借助于變分原理或加權余量法，將微分方程離散求解[3]。采用不同的權函數和插值函數形式，便構成不同的有限元方法。

2.3 有限元法的計算格式

根據所采用的權函數和插值函數的不同，有限元方法也分為多種計算格式。從權函數的選擇來說，有配置法、矩量法、最小二乘法和伽遼金法；從計算單元網格的形狀來劃分，有三角形網格、四邊形網格和多邊形網格；從插值函數的精度來劃分，又分為線性插值函數和高次插值函數等。不同的組合同樣構成不同的有限元計算格式。插值函數一般由不同次冪的多項式組成，但也有采用三角函數或指數函數組成的乘積表示，但最常用的多項式插值函數。有限元插值函數分為兩大類，一類只要求插值多項式本身在插值點取已知值，稱為拉格朗日（Lagrange）多項式插值；另一種不僅要求插值多項式本身，還要求它的導數值在插值點取已知值，稱為哈密特（Hermite）多項式插值。在二維有限元中，三角形單元應用的最早，近來四邊形等參元的應用也越來越廣[4]。對于二維三角形和四邊形電源單元，常采用的插值函數為有Lagrange插值直角坐標系中的線性插值函數及二階或更高階插值函數、面積坐標系中的線性插值函數、二階或更高階插值函數等。

2.4 有限元法的解題步驟

對于有限元方法，其基本思路和解題步驟可歸納為：

（1）建立積分方程：根據變分原理或方程余量與權函數正交化原理，建立與微分方程初邊值問題等價的積分表達式，這是有限元法的出發點。

（2）區域單元剖分：根據求解區域的形狀及實際問題的物理特點，將區域剖分為若干相互連接、不重疊的單元[5]。區域單元劃分是采用有限元方法的前期準備工作，這部分工作量比較大，除了給計算單元和節點進行編號和確定相互之間的關系之外，還要表示節點的位置坐標，同時還需要列出自然邊界和本質邊界的節點序號和相應的邊界值。

（3）確定單元基函數：根據單元中節點數目及對近似解精度的要求，選擇滿足一定插值條件的插值函數作為單元基函數。有限元方法中的基函數是在單元中選取的，由于各單元具有規則的幾何形狀，在選取基函數時可遵循一定的法則。

（4）單元分析：將各個單元中的求解函數用單元基函數的線性組合表達式進行逼近；再將近似函數代入積分方程，并對單元區域進行積分，可獲得含有待定系數（即單元中各節點的參數值）的代數方程組，稱為單元有限元方程。

（5）總體合成：在得出單元有限元方程之后，將區域中所有單元有限元方程按一定法則進行累加，形成總體有限元方程。

（6）邊界條件的處理：一般邊界條件有三種形式，分為本質邊界條件（狄里克雷邊界條件）、自然邊界條件（黎曼邊界條件）、混合邊界條件（柯西邊界條件）。對于自然邊界條件，一般在積分表達式中可自動得到滿足。對于本質邊界條件和混合邊界條件，需按一定法則對總體有限元方程進行修正滿足。

（7）解有限元方程：根據邊界條件修正的總體有限元方程組，是含所有待定未知量的封閉方程組，采用適當的數值計算方法求解，可求得各節點的函數值[6]。

3 有限元法在MATLAB中的實現（Reality of finite

element numerical method in MATLAB）

在用有限元法時，MATLAB編程的方法與其他語言如FORTRAN、C語言相似，但用MATLAB編程會更簡單一些。以下詳細介紹了用有限元方法計算十字形截面的方形同軸電纜線內的電勢，該同軸電纜線內的電勢值在每一點都不一樣，無法用解析解來描述，所以用了有限元的方法解數值解，最后畫出了等勢線，對程序稍加修改，也可以用表面圖反映數值解的分布。

（1）問題描述

有橫截面如圖1所示的電纜芯線，電纜的外導體的電勢為零，芯線的電勢為10V。求該同軸電纜線內的電勢[7]。

（2）定解問題

這可以看成是二維問題，由于問題的對稱性，只要解第一象限（四分之一區域）就可以，根據問題建立如下的方程組：

（3）區域單元剖分

根據本問題的性質，將區域劃分為三角形，網格的劃分如圖2所示。

4 結論（Conclusion）

科學計算在各門自然科學（物理學、氣象學、地質學和生命科學等）和技術科學與工程科學（核技術、石油勘探、航空航天和大型土木工程等）中起著越來越重要的作用，在很多重要領域中成為不可缺少的工具。而科學與工程計算中最重要的內容就是求解在科學研究和工程技術中出現的各種各樣的偏微分方程或方程組。解偏微分方程已經成為科學與工程計算的核心內容，包括一些大型的計算和很多已經成為常規的計算。原則上，可以用FORTRAN或C語言來完成這些計算，但很少有人這樣去做，原因是成本太高，編程太復雜。而MATLAB是一種用于算法開發、數據可視化、數據分析以及數值計算的高級技術計算語言和交互式環境[8]。所以使用MATLAB，可以較使用傳統的編程語言（如C、C++和Fortran）更快地解決技術計算問題。

參考文獻（References）

[1] 陸金甫，關治.偏微分方程數值解法（第2版）[M].北京：清華大

學出版社，2003.

[2] 張義寬，張曉濱.計算機圖形學[M].西安：西安電子科技大學出

版社，2004.

[3] 周俊明，林群.高次三角形有限元外推的探討[J].數學實踐與認

識，2008，38（5）：99-106.

[4] 周俊明，林群.高次三角形有限元外推的進一步研究[J].數學

的實踐與認識，2008，38（16）：192-197.

[5] 李長河，等.基于網格數據圖的自適應細分及邊緣提取[J].計

算機工程與應用，2003，39（24）：94-96；109.

[6] 楊曉松，顧元憲.有限元網格體繪制中的剖切算法[J].中國圖

象圖形學報，2002，（1）：55-62.

[7] 彭芳麟.數物理方程的MATLAB解法與可視化[M].北京：清華

大學出版社，2004.

[8] 于萬波.基于MATLAB的計算機圖形與動畫技術[M].北京：清

華大學出版社，2007.

作者簡介：

馮桂蓮（1979-），女，碩士，副教授.研究領域：計算機專業

教學，計算機軟件與理論，計算機圖形學.