有限元數(shù)值解法在MATLAB中的實(shí)現(xiàn)及可視化

摘 要：偏微分方程的數(shù)值解法在數(shù)值分析中占有很重要的地位，很多科學(xué)技術(shù)問(wèn)題的數(shù)值計(jì)算包括了偏微分方程的數(shù)值解問(wèn)題。在學(xué)習(xí)初等函數(shù)時(shí)，總是先畫(huà)出它們的圖形，因?yàn)閳D形能幫助我們了解函數(shù)的性質(zhì)。而對(duì)于偏微分方程，畫(huà)出它們的圖形并不容易，尤其是沒(méi)有解析解的偏微分方程，畫(huà)圖就顯得更加不容易了。為了從偏微分方程的數(shù)學(xué)表達(dá)式中看出其所表達(dá)的圖形、函數(shù)值與自變量之間的關(guān)系，通過(guò)MATLAB編程，用有限元數(shù)值解法求解了偏微分方程，并將其結(jié)果可視化。

關(guān)鍵詞：偏微分方程；MATLAB；有限元法；可視化

中圖分類(lèi)號(hào)：TP311.12 文獻(xiàn)標(biāo)識(shí)碼：A

1 引言（Introduction）

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

在學(xué)習(xí)初等函數(shù)時(shí)，總是先畫(huà)出它們的圖形，因?yàn)閳D形能幫助我們了解函數(shù)的性質(zhì)。而對(duì)于偏微分方程，畫(huà)出它們的圖形并不容易，尤其是沒(méi)有解析解的偏微分方程，畫(huà)圖就顯得更加不容易了。所以本文主要研究如何用MATLAB數(shù)值求解偏微分方程，并將其數(shù)值解繪制成三維圖形的形式，從而可以從復(fù)雜的數(shù)學(xué)表達(dá)式中看出其所表達(dá)的圖像、函數(shù)值與自變量之間的關(guān)系[2]。

2 有限元法（Finite element method）

2.1 有限元法概述

有限元法的基本思想是將結(jié)構(gòu)離散化，用有限個(gè)容易分析的單元來(lái)表示復(fù)雜的對(duì)象，單元之間通過(guò)有限個(gè)節(jié)點(diǎn)相互連接，然后根據(jù)變形協(xié)調(diào)條件綜合求解。由于單元的數(shù)目是有限的，節(jié)點(diǎn)的數(shù)目也是有限的，所以稱(chēng)為有限元法。

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

2.2 有限元法的基本思想

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

有限元方法的基礎(chǔ)是變分原理和加權(quán)余量法，其基本求解思想是把計(jì)算域劃分為有限個(gè)互不重疊的單元，在每個(gè)單元內(nèi)，選擇一些合適的節(jié)點(diǎn)作為求解函數(shù)的插值點(diǎn)，將微分方程中的變量改寫(xiě)成由各變量或其導(dǎo)數(shù)的節(jié)點(diǎn)值與所選用的插值函數(shù)組成的線(xiàn)性表達(dá)式，借助于變分原理或加權(quán)余量法，將微分方程離散求解[3]。采用不同的權(quán)函數(shù)和插值函數(shù)形式，便構(gòu)成不同的有限元方法。

2.3 有限元法的計(jì)算格式

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

2.4 有限元法的解題步驟

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

（1）建立積分方程：根據(jù)變分原理或方程余量與權(quán)函數(shù)正交化原理，建立與微分方程初邊值問(wèn)題等價(jià)的積分表達(dá)式，這是有限元法的出發(fā)點(diǎn)。

（2）區(qū)域單元剖分：根據(jù)求解區(qū)域的形狀及實(shí)際問(wèn)題的物理特點(diǎn)，將區(qū)域剖分為若干相互連接、不重疊的單元[5]。區(qū)域單元?jiǎng)澐质遣捎糜邢拊椒ǖ那捌跍?zhǔn)備工作，這部分工作量比較大，除了給計(jì)算單元和節(jié)點(diǎn)進(jìn)行編號(hào)和確定相互之間的關(guān)系之外，還要表示節(jié)點(diǎn)的位置坐標(biāo)，同時(shí)還需要列出自然邊界和本質(zhì)邊界的節(jié)點(diǎn)序號(hào)和相應(yīng)的邊界值。

（3）確定單元基函數(shù)：根據(jù)單元中節(jié)點(diǎn)數(shù)目及對(duì)近似解精度的要求，選擇滿(mǎn)足一定插值條件的插值函數(shù)作為單元基函數(shù)。有限元方法中的基函數(shù)是在單元中選取的，由于各單元具有規(guī)則的幾何形狀，在選取基函數(shù)時(shí)可遵循一定的法則。

（4）單元分析：將各個(gè)單元中的求解函數(shù)用單元基函數(shù)的線(xiàn)性組合表達(dá)式進(jìn)行逼近；再將近似函數(shù)代入積分方程，并對(duì)單元區(qū)域進(jìn)行積分，可獲得含有待定系數(shù)（即單元中各節(jié)點(diǎn)的參數(shù)值）的代數(shù)方程組，稱(chēng)為單元有限元方程。

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

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

（7）解有限元方程：根據(jù)邊界條件修正的總體有限元方程組，是含所有待定未知量的封閉方程組，采用適當(dāng)?shù)臄?shù)值計(jì)算方法求解，可求得各節(jié)點(diǎn)的函數(shù)值[6]。

3 有限元法在MATLAB中的實(shí)現(xiàn)（Reality of finite

element numerical method in MATLAB）

在用有限元法時(shí)，MATLAB編程的方法與其他語(yǔ)言如FORTRAN、C語(yǔ)言相似，但用MATLAB編程會(huì)更簡(jiǎn)單一些。以下詳細(xì)介紹了用有限元方法計(jì)算十字形截面的方形同軸電纜線(xiàn)內(nèi)的電勢(shì)，該同軸電纜線(xiàn)內(nèi)的電勢(shì)值在每一點(diǎn)都不一樣，無(wú)法用解析解來(lái)描述，所以用了有限元的方法解數(shù)值解，最后畫(huà)出了等勢(shì)線(xiàn)，對(duì)程序稍加修改，也可以用表面圖反映數(shù)值解的分布。

（1）問(wèn)題描述

有橫截面如圖1所示的電纜芯線(xiàn)，電纜的外導(dǎo)體的電勢(shì)為零，芯線(xiàn)的電勢(shì)為10V。求該同軸電纜線(xiàn)內(nèi)的電勢(shì)[7]。

（2）定解問(wèn)題

這可以看成是二維問(wèn)題，由于問(wèn)題的對(duì)稱(chēng)性，只要解第一象限（四分之一區(qū)域）就可以，根據(jù)問(wèn)題建立如下的方程組：

（3）區(qū)域單元剖分

根據(jù)本問(wèn)題的性質(zhì)，將區(qū)域劃分為三角形，網(wǎng)格的劃分如圖2所示。

4 結(jié)論（Conclusion）

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

參考文獻(xiàn)（References）

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

學(xué)出版社，2003.

[2] 張義寬，張曉濱.計(jì)算機(jī)圖形學(xué)[M].西安：西安電子科技大學(xué)出

版社，2004.

[3] 周俊明，林群.高次三角形有限元外推的探討[J].數(shù)學(xué)實(shí)踐與認(rèn)

識(shí)，2008，38（5）：99-106.

[4] 周俊明，林群.高次三角形有限元外推的進(jìn)一步研究[J].數(shù)學(xué)

的實(shí)踐與認(rèn)識(shí)，2008，38（16）：192-197.

[5] 李長(zhǎng)河，等.基于網(wǎng)格數(shù)據(jù)圖的自適應(yīng)細(xì)分及邊緣提取[J].計(jì)

算機(jī)工程與應(yīng)用，2003，39（24）：94-96；109.

[6] 楊曉松，顧元憲.有限元網(wǎng)格體繪制中的剖切算法[J].中國(guó)圖

象圖形學(xué)報(bào)，2002，（1）：55-62.

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

大學(xué)出版社，2004.

[8] 于萬(wàn)波.基于MATLAB的計(jì)算機(jī)圖形與動(dòng)畫(huà)技術(shù)[M].北京：清

華大學(xué)出版社，2007.

作者簡(jiǎn)介：

馮桂蓮（1979-），女，碩士，副教授.研究領(lǐng)域：計(jì)算機(jī)專(zhuān)業(yè)

教學(xué)，計(jì)算機(jī)軟件與理論，計(jì)算機(jī)圖形學(xué).