黃鵬展

（新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,烏魯木齊 830046）

眾所周知,大量實(shí)際中的工程問(wèn)題都可以用偏微分方程來(lái)描述,但有許多偏微分方程定解問(wèn)題的解不能以實(shí)用的解析形式來(lái)表出,從而無(wú)法得到這些方程的準(zhǔn)確解,以定量地描述客觀過(guò)程。有限元課程是數(shù)值分析的后繼課程,是計(jì)算數(shù)學(xué)專業(yè)學(xué)生必修的專業(yè)課程。課程中的有限元方法在數(shù)值計(jì)算的理論基礎(chǔ)上,構(gòu)造求解的近似逼近方法。該課程通過(guò)一些典型、常用、有效的數(shù)值方法來(lái)解釋有限元方法的基本思想,使學(xué)生了解如何在計(jì)算機(jī)上應(yīng)用這些數(shù)值方法求解一個(gè)偏微分方程定解問(wèn)題。同時(shí),通過(guò)對(duì)有限元基本概念及基本理論的學(xué)習(xí),使學(xué)生的理論分析能力得到一定的訓(xùn)練,并通過(guò)學(xué)習(xí)與實(shí)踐使得學(xué)生熟練掌握及運(yùn)用有限元方法來(lái)數(shù)值求解偏微分方程。

而有限元實(shí)驗(yàn)課程的設(shè)立,則能更好地為學(xué)生學(xué)習(xí)有限元服務(wù),使得學(xué)生能夠更好、更深刻地了解有限元方法,體會(huì)有限元方法在數(shù)值求解偏微分方程時(shí)的魅力。本文主要結(jié)合筆者的親身教學(xué)經(jīng)歷,從有限元的Matlab工具包、有限元軟件Free?fem++以及Matlab編程等三個(gè)方面來(lái)探索有限元實(shí)驗(yàn)課程的教學(xué)方法。這三個(gè)方面從簡(jiǎn)到難、從特定偏微分方程類型求解到較一般偏微分方程求解,循序漸進(jìn)地引導(dǎo)學(xué)生進(jìn)行有限元方法的實(shí)踐和數(shù)值的實(shí)驗(yàn)。

1 有限元方法

本文以文獻(xiàn)〔1〕中的一個(gè)偏微分方程為例,簡(jiǎn)單地介紹基本的有限元方法。

考慮二維區(qū)域上的Poisson方程〔2〕第三邊值問(wèn)題

這里的Ω∈R2是具有光滑邊界的有界凸區(qū)域,n是?Ω上的外法線方向,α和g是?Ω上的連續(xù)函數(shù),α=α(x,y)≥0,△是拉普拉斯算子,f是作用力。為簡(jiǎn)便起見,在后面的計(jì)算部分,令 f=g=α=1,計(jì)算區(qū)域Ω取成單位正方形。

我們知道有限元方法基于變分原理,故先給出方程（1）的變分形式。即第一步,

第二步,給出計(jì)算區(qū)域的網(wǎng)格剖分。設(shè)h是一個(gè)趨近于0的正參數(shù),記網(wǎng)格剖分為Kh,它由一些三角單元K組成。

第三步,試探函數(shù)空間的構(gòu)造。即選取有限元子空間

其中, φi（i=1,2,…,N）是基函數(shù),滿足

求使得，D(uh,vh)=F(vh),?vh∈ Xh，

其中ui=u(xi,yi)。

第四步,有限元離散。生成單元?jiǎng)偠染仃嚭蛦卧奢d向量,通過(guò)單元上局部編碼與整體編碼的對(duì)應(yīng)關(guān)系得到的擴(kuò)張矩陣,生成總剛度矩陣和總荷載向量,然后對(duì)約束條件進(jìn)行處理,最后得到有限元方程組。由方程組的求解方法解出ui,結(jié)合基函數(shù)φi得到有限元近似解uh。

2 有限元上機(jī)實(shí)驗(yàn)

有限元課程是一門計(jì)算數(shù)學(xué)專業(yè)中較難的課程,故而,我們首先對(duì)問(wèn)題（1）使用較為簡(jiǎn)單Matlab工具包來(lái)實(shí)現(xiàn)數(shù)值模擬,鑒于它的可視化操作,學(xué)生易于上手。從而,增強(qiáng)他們對(duì)這門較難課程學(xué)習(xí)的信心,克服學(xué)習(xí)前的恐懼心理和由此產(chǎn)生的怠學(xué)狀態(tài)。通過(guò)Matlab工具包對(duì)特定方程的模擬,能夠使學(xué)生了解所求得方程的解到底是什么樣子,初步了解有限元解法的過(guò)程。接著,采用有限元軟件Freefem++來(lái)數(shù)值求解問(wèn)題（1）,該有限元軟件只需要寫出有限元的變分形式即可。最后,按照有限元編程的框架圖,使用Matlab來(lái)編制有限元求解（1）的程序。

2.1 Matlab工具包 根據(jù)文獻(xiàn)〔3〕中的講解,筆者以（1）為例,對(duì)Matlab工具包求解偏微分方程進(jìn)行介紹。

第一步,打開Matlab軟件,在命令框中輸入pdetool,便可啟動(dòng)一個(gè)可視化界面。接著,完成求解區(qū)域的構(gòu)造。

第二步,選取邊界,打開Boundary Condition對(duì)話框,輸入（1）中的第三類邊界條件。

第三步,選擇PDE Mode命令,設(shè)置方程的類型,這里為橢圓型方程,接著輸入（1）中偏微分方程的系數(shù)。

第四步,利用Initialize Mesh命令,對(duì)第一步構(gòu)造的求解區(qū)域進(jìn)行網(wǎng)格剖分。如果網(wǎng)格不夠密,可以利用Refine Mesh命令進(jìn)行加密。

第五步,選擇Solve PDE命令,求解設(shè)定的偏微分方程適定問(wèn)題,給出數(shù)值解uh的圖形,見圖1。

圖1 用Matlab工具pdetool得到（1）的有限元解

Matlab工具包可以對(duì)空間二維問(wèn)題快速高效求解。由于它的界面可視化,學(xué)生只需要使用用戶界面,構(gòu)造出求解區(qū)域,輸入方程的類型和邊界條件,自動(dòng)生成網(wǎng)格,不需要人工設(shè)定基函數(shù),就可以顯示出數(shù)值解的結(jié)果。故而,它是一種利用有限元方法數(shù)值求解偏微分方程的簡(jiǎn)單有效工具。但是這個(gè)可視化界面只是對(duì)特定的一些方程,并且邊界條件的設(shè)置也很固定。最主要的是不能夠自己選擇基函數(shù)或者構(gòu)造基函數(shù),而有限元的魅力主要正是在于它所擁有的各式各樣的基函數(shù),并可以根據(jù)自己的需要進(jìn)行構(gòu)造。

故而,Matlab工具包只適合有限元方法的初學(xué)者,使他們較好較快地了解求解偏微分方程的有限元方法。在這方面,自然地它有不可替代的優(yōu)勢(shì)。因此,筆者在有限元實(shí)驗(yàn)課程教學(xué)時(shí),最初幾節(jié)課,主要就是給學(xué)生們講解Matlab工具包求解偏微分方程。

2.2 有限元軟件 FreeFem++ FreeFem++〔4〕是一款數(shù)值求解偏微分方程的有限元軟件。就像它名字所說(shuō)的,它是一款免費(fèi)的軟件,可以直接在http://www.freefem.org這個(gè)網(wǎng)站上下載安裝包,并自帶它的教程。FreeFem++不是一個(gè)工具包,而是一個(gè)完整的有限元產(chǎn)品,具有內(nèi)存需求少、計(jì)算速度快、易學(xué)好用、功能強(qiáng)大、應(yīng)用面廣等特點(diǎn)〔5-8〕。故而,在有限元實(shí)驗(yàn)課程隨后的課時(shí)中,筆者將講授這部分的內(nèi)容。

仍以（1）為例,對(duì)有限元軟件FreeFem++求解偏微分方程進(jìn)行介紹。以下部分就是用該軟件實(shí)現(xiàn)數(shù)值求解的FreeFem++程序,附程序注釋。

real n=80;//設(shè)置網(wǎng)格尺度

mesh Th=square（n,n）;//生成網(wǎng)格

fespace Xh（Th,P2）;//定義有限元空間

Xh u,v;//聲明試探函數(shù)和檢驗(yàn)函數(shù)

func f=1;//選取作用力函數(shù)

solve Poisson（u,v）=//定義偏微分方程,即變分形式

·馬改戶 曹春生 何鄂 陳云崗 黎明 郭北平楊參軍 白羽平 冷軍 朱春林等10位雕塑、油畫大家力挺學(xué)術(shù)邀請(qǐng)展

int2d（Th）（dx（u）*dx（v）+dy（u）*dy（v））//雙線性項(xiàng)形式

+int1d（Th）（u*v）//D（u,v）的第二部分,系數(shù)為1

-int2d（Th）（v）//右端項(xiàng)

-int1d（Th）（v）;//F（v）的第二部分,g取1

plot（u,ps=“fig2.eps”）;//計(jì)算結(jié)果的可視化這樣就完成了采用FreeFem++軟件求解（1）的全部過(guò)程,數(shù)值結(jié)果見圖2。

圖2 利用FreeFem++軟件得到（1）的有限元解

由以上的程序我們可以看到,使用該軟件進(jìn)行編程數(shù)值求解偏微分方程時(shí),只需要寫出原問(wèn)題的變分問(wèn)題即可,無(wú)需進(jìn)行單元局部和總體的編號(hào)、單元?jiǎng)偠染仃嚭蛦卧奢d向量的計(jì)算、總剛度矩陣和總荷載向量的組成、有限元線性方程組的求解等具體過(guò)程。更加主要的是采用“fespace Xh（Th,P2）”這段語(yǔ)句就完成了有限元空間的定義,并且使用“int2d（Th）”這個(gè)命令就實(shí)現(xiàn)了數(shù)值積分的功能。寥寥數(shù)行即完成了對(duì)Poisson方程第三類邊值問(wèn)題的有限元方法求解。對(duì)比使用C或者M(jìn)atlab編程需幾百行的代碼,FreeFem++程序大大節(jié)約了編程時(shí)間。

但是,由于它是一款免費(fèi)軟件,其底層語(yǔ)言沒(méi)有公開,繼而我們無(wú)法了解最本質(zhì)的內(nèi)容。故而,它有一定的局限性,無(wú)法對(duì)底層的代碼進(jìn)行修改,二次編程。并且,只能選取它所給定的一些基函數(shù)和有限元空間。因此,它適合一些想使用有限元方法做一些初步研究的研究生。如果想要驗(yàn)證更加原創(chuàng)的一些有限元算法,筆者認(rèn)為應(yīng)該利用C或Matlab或者Fortran語(yǔ)言來(lái)實(shí)現(xiàn)。

2.3基于Matlab的有限元編程 鑒于以上考慮,筆者在這門課程的剩余課時(shí)里,重點(diǎn)講解基于Matlab的有限元方法編程。主要是根據(jù)文獻(xiàn)〔9-10〕中的程序框圖進(jìn)行程序編制,仍舊以（1）為例。

M文件1:構(gòu)造基函數(shù) function J=Base（m,n），m,n分別是x,y軸等分點(diǎn)數(shù)。

M文件2:構(gòu)造矩陣function A=juzhen（n,m），即構(gòu)造求解未知節(jié)點(diǎn)所滿足的矩陣。這里為重點(diǎn)部分,矩陣由總剛度矩陣并利用數(shù)值積分得到。

M文件3:當(dāng)f=1時(shí)矩陣的右端項(xiàng)function b=youduanxiang（n,m）。

M文件4:將求出的矩陣變換為一般的矩陣function B=bh（n,m），即將矩陣A轉(zhuǎn)換成與理論分析時(shí)計(jì)算的矩陣一樣。

M文件5:右端項(xiàng)為1時(shí),用CG方法求解矩陣function u=gonge（n,m,detla,max），這里,max為最大循環(huán)次數(shù),detla=1.0E-6為容許誤差。

M文件6:右端項(xiàng)為1時(shí),用有限元方法得到的解 U=qiujie（n,m,detla,max）。

再利用一個(gè)主文件調(diào)用以上的M文件,給定所設(shè)的參數(shù),就可以得到（1）的有限元方法數(shù)值解,見圖3。

圖3 利用Matlab編程得到（1）的有限元解

故而,利用文獻(xiàn)〔9〕、〔10〕里面的程序框架圖,我們可以編制有限元的Matlab程序,從而得到想要的數(shù)值結(jié)果。直接編程的主要優(yōu)點(diǎn)就是算法驗(yàn)證的靈活性,可以任意構(gòu)造希望的基函數(shù)及有限元空間。但是它所花費(fèi)的時(shí)間最長(zhǎng),因編制程序以及調(diào)試程序等都得消耗大量的時(shí)間。

3 結(jié)論

本文結(jié)合筆者的教學(xué)經(jīng)歷,對(duì)有限元實(shí)驗(yàn)課程教學(xué)進(jìn)行了探索。針對(duì)數(shù)值求解偏微分方程,主要關(guān)于Matlab工具包、有限元軟件Freefem++以及Matlab有限元編程這三個(gè)方面展開討論。具體介紹了這三種方法的使用和編制,給出了每種方法的優(yōu)缺點(diǎn)。

總之,對(duì)有限元實(shí)驗(yàn)課程的教學(xué)要循序漸進(jìn)。Matlab工具包最為簡(jiǎn)單適合初學(xué)者使用,但是對(duì)于求解的問(wèn)題具有局限性;有限元軟件Freefem++次之,但是對(duì)基函數(shù)選取與構(gòu)造等方面的靈活性欠缺,故它適合成為研究生用來(lái)做一些初步研究的有效工具;Matlab有限元編程相比之下最為困難與耗時(shí),但它是驗(yàn)證一些有限元原創(chuàng)算法的必備工具。

〔1〕陸金甫,關(guān)治.偏微分方程數(shù)值解〔M〕.2版.北京:清華大學(xué)出版社,2004.

〔2〕劉相國(guó),謝如龍,郝江鋒.二維泊松方程的交替方向迭代法〔J〕.大理學(xué)院學(xué)報(bào),2010,9（10）:1-5.

〔3〕陸君安,尚濤,謝進(jìn),等.偏微分方程的MATLAB解法〔M〕.武漢:武漢大學(xué)出版社,2001.

〔4〕Hecht F.FreeFem++〔EB/OL〕.（2012-03-16）〔2012-10-10〕.http://www.freefem.org/ff++.

〔5〕尚月強(qiáng).FreeFem++與偏微分方程有限元數(shù)值計(jì)算〔J〕.貴州師范學(xué)院學(xué)報(bào),2012,27（12）:1-5.

〔6〕尚月強(qiáng).基于MPI+FreeFem++的有限元并行計(jì)算〔J〕.軟件導(dǎo)刊,2012,11（2）:27-29.

〔7〕杜世森.基于FreeFEM++的一階有限元解法〔D〕.長(zhǎng)沙:湖南大學(xué),2011.

〔8〕詹陽(yáng)烈,莊春剛，熊振華，等.基于水平集方法與Free?FEM的彈性結(jié)構(gòu)拓?fù)鋬?yōu)化〔J〕.中國(guó)機(jī)械工程,2009,20（11）:1326-1331.

〔9〕李開泰,黃艾香,黃慶懷.有限元方法及其應(yīng)用〔M〕.北京:科學(xué)出版社,2006.

〔10〕馬逸塵,梅立泉,王阿霞.偏微分方程現(xiàn)代數(shù)值方法〔M〕.北京:科學(xué)出版社,2006.