橢圓型偏微分方程的三角形單元有限元的數(shù)值解法

俞 輝， 楊臻豪， 奉繼明， 瞿少成

(1.三峽大學(xué) 理學(xué)院， 湖北 宜昌 443000; 2.華中師范大學(xué) 電子信息工程系， 武漢 430079)

?

橢圓型偏微分方程的三角形單元有限元的數(shù)值解法

俞 輝1*， 楊臻豪1， 奉繼明1， 瞿少成2

(1.三峽大學(xué) 理學(xué)院， 湖北 宜昌 443000; 2.華中師范大學(xué) 電子信息工程系， 武漢 430079)

本文以橢圓型偏微分方程問題為背景研究三角形單元的有限元方法.隨著計(jì)算機(jī)圖形學(xué)的發(fā)展，三角形單元?jiǎng)澐秩〉昧司薮蟮某删停梢缘玫劫|(zhì)量非常好的三角形單元，進(jìn)而提高偏微分方程的有限元方法的數(shù)值解的精度.本文采用節(jié)點(diǎn)增量算法，對(duì)問題區(qū)域進(jìn)行三角形單元?jiǎng)澐郑玫降娜切螁卧獫M足Delaunay條件，再對(duì)三角形單元的所有節(jié)點(diǎn)采用自適應(yīng)編號(hào)，最后運(yùn)用三角形單元的有限元方法得到橢圓型偏微分方程的數(shù)值解.通過數(shù)值實(shí)驗(yàn)，得出相比傳統(tǒng)的三角形單元的有限元方法，本文的三角形單元的有限元方法減小了舍入誤差，提高了計(jì)算精度.

Delaunay三角形單元； 自適應(yīng)編號(hào)； 舍入誤差

隨著經(jīng)濟(jì)的發(fā)展，工程的各個(gè)領(lǐng)域都呈現(xiàn)蓬勃發(fā)展的趨勢(shì).但是在發(fā)展的同時(shí)遇到了很多的工程問題.很大一部分的工程問題通過數(shù)學(xué)建模[1]，可以轉(zhuǎn)化成偏微分方程的問題.通過對(duì)流體力學(xué)的分析，應(yīng)用數(shù)學(xué)知識(shí)，可以將流體力學(xué)中的問題轉(zhuǎn)化成偏微分方程，即Navier-Stokes[2]方程，通過方程的分析和求解其數(shù)值解，足以描述邊界層、渦流等多種現(xiàn)象.在基于股票的衍生證券市場(chǎng)上，對(duì)買入期權(quán)、賣出期權(quán)、股票的買入價(jià)格、賣出價(jià)格的分析，布萊爾(Black[3])和休爾斯(Scholes[3])得出了歐式期權(quán)的無套利價(jià)格的倒向偏微分方程.可以看出偏微分方程在很多領(lǐng)域都涉及到了.

有限元方法[4]是求解偏微分方程數(shù)值解的重要方法.有限元求解偏微分方程時(shí)，一個(gè)重要的步驟是[5]，運(yùn)用數(shù)學(xué)知識(shí)將偏微分方程轉(zhuǎn)化為積分方程.對(duì)偏微分方程的原問題區(qū)域進(jìn)行單元?jiǎng)澐諿6-8]，如四邊形單元?jiǎng)澐趾腿菃卧獎(jiǎng)澐?20世紀(jì)80年代三角網(wǎng)格剖分在眾多領(lǐng)域引起了廣泛的關(guān)注.英國Bath大學(xué)數(shù)學(xué)分校的P.J Green和R. Sibson等從數(shù)學(xué)的角度對(duì)三角網(wǎng)格剖分進(jìn)行了研究[9].Bath大學(xué)數(shù)學(xué)分校的A. Bowyer和澳大利亞悉尼大學(xué)Geology and Geophyics系的D.F.Watson于1981年發(fā)表了文章[10]，提出各不相同的三角網(wǎng)格構(gòu)造方法，但最終得到的結(jié)果都是以Dirichlet/Voronoi圖為理論基礎(chǔ)的Delaunay三角網(wǎng)格.

本文將先進(jìn)的三角形單元網(wǎng)格劃分運(yùn)用到有限元方法中.對(duì)問題區(qū)域采用成熟的節(jié)點(diǎn)增量算法[11]進(jìn)行單元剖分，得到的三角形單元都是滿足Delaunay條件的三角形.在運(yùn)用有限元方法時(shí)，單元的“剖分”、“設(shè)置”、“編號(hào)”將影響有限元方程系數(shù)矩陣的形狀，最后影響結(jié)果的計(jì)算精度.本文對(duì)問題區(qū)域進(jìn)行三角形單元剖分后，對(duì)單元中的所有節(jié)點(diǎn)進(jìn)行的編號(hào)，計(jì)算單元上系數(shù)矩陣和組裝單元矩陣時(shí)，對(duì)單元的節(jié)點(diǎn)采用自適應(yīng)編號(hào)排序，得到單元節(jié)點(diǎn)關(guān)系.最后從數(shù)值仿真例子可以得出，本文的有限元方法，比傳統(tǒng)的有限元方法，計(jì)算精度更高.

1 有限元構(gòu)思

考慮以下經(jīng)典的二階橢圓型方程的邊值問題[12]：

(1)

其中，Ω是平面有界區(qū)域， Γ是區(qū)域的邊界.

.

引理[5]設(shè)F是一個(gè)微分算子，則

(2)

通過上述引理(2)，可以將經(jīng)典的橢圓型偏微分方程轉(zhuǎn)化成積分方程了，即在邊值問題中的偏微分方程兩邊分別同時(shí)乘以v，再對(duì)兩邊進(jìn)行積分，從而將偏微分方程問題，轉(zhuǎn)化如下成積分方程：

(3)

其中

從而將式(3)化簡得到下面式子

(4)

記

于是得到如下式子

(5)

考慮到二階橢圓型偏微分方程(1)的邊界條件，即

其中，Γ是問題域的邊界.所以可得到

于是橢圓型偏微分方程邊值問題可以轉(zhuǎn)化為

(6)

上式(6)稱為橢圓型偏微分方程邊值問題的積分形式.

從問題的積分形式可以看出，二階橢圓型偏微分方程的條件是u的二階偏導(dǎo)存在，通過運(yùn)用Green格林公式，將u的二階偏導(dǎo)問題轉(zhuǎn)化為一階偏導(dǎo)問題，這就使得原問題的解u的光滑性的要求大為下降，只要u和偏導(dǎo)u′x，u′y存在并且可積就可以了.將有限元方法應(yīng)用于二階橢圓型偏微分方程的邊值問題，其本質(zhì)是將有限元方法應(yīng)用于積分形式的方程，故本文以后將放棄二階橢圓型偏微分方程的邊值問題的原形式，而直接從積分形式出發(fā)來研究.

設(shè)Pm為次數(shù)不超過m的多項(xiàng)式空間，則本文所要用到的多項(xiàng)式樣條函數(shù)空間，可以寫為

(7)

(8)

2 三角形單元的有限元方法實(shí)現(xiàn)

有限元方法中一個(gè)重要的思想是，將問題區(qū)域進(jìn)行劃分，得到一個(gè)一個(gè)單元區(qū)域，對(duì)單元區(qū)域應(yīng)用有限元方法得到單元矩陣，通過單元區(qū)域之間的關(guān)系，組裝單元矩陣得到線性方程組的系數(shù)矩陣.故三角形單元?jiǎng)澐諿14]包括2大部分，一是對(duì)問題區(qū)域進(jìn)行三角形劃分；二是建立單元與單元之間的關(guān)系即單元關(guān)系矩陣.

2.1三角形單元?jiǎng)澐?/p>

本文仿真例子的問題區(qū)域是Ω={(x，y)|x2+y2=1}.取步長hx=hy=0.5，對(duì)圓形區(qū)域進(jìn)行三角形單元?jiǎng)澐?在三角形單元?jiǎng)澐謺r(shí)，要避免畸形單元，即三角單元中存在角度非常小的單元.為了避免畸形單元，本文采用成熟的節(jié)點(diǎn)增量算法，使三角形單元滿足Delaunay三角化特性[11].得圖1.

圖1 圓形區(qū)域三角形單元?jiǎng)澐諪ig.1 Circular area triangle unit division

2.2單元變換

圖2 單元變換Fig.2 Unit transform

設(shè)

2.3單元插值函數(shù)

設(shè)

運(yùn)用Matlab編程計(jì)算，可得：c1=2，c2=0，c3=0，c4=-1，c5=0，c6=0

Ni=2ξ2-ξ.

同理可得

Nj=2η2-η，

Nk=2ξ2+2η2+4ξη-3ξ-3η+1，

Nr=4ξη，

Nq=4ξ-4ξη-4ξ2， Np=4η-4ξη-4η2.

至此就建立了標(biāo)準(zhǔn)三角單元上Lagrange二次型的6個(gè)基函數(shù)了，故Lagrange二次型插值函數(shù)如下

Npup+Nquq+Nrur.

2.4有限元方法的實(shí)現(xiàn)

在建立有限元方程[12]時(shí)，用的是Lagrange二次型插值函數(shù)，直接從微分方程的積分方程出發(fā)，積分方程的形式如下

因?yàn)長agrange二次型插值函數(shù)為

Npup+Nquq+Nrur,

Npvp+Nqvq+Nrvr,

將插值函數(shù)代入得

aijuivj+aikuivk+aipuivp+…+

ariurvi+arjurvj+arkurvk+arpurvp+

μijuivj+μiruivr+…+μriurvi+

djvj+dkvk+dpvp+dqvq+drvr],

其中,γij是邊界Γh的一條邊，長度為lij，并且在γij上，ξ+η=1.

q·NmNl]Jdξdη,m，l=i，j，k，p，q，r,

令

故

上式即為三角單元的有限元方程，通過解線性方程就可以求解出橢圓型偏微分方程有限元的數(shù)值解了.

2.5單元自適應(yīng)編號(hào)排序

下圖3和圖4是2種單元變換.

圖3 單元變換1Fig.3 1 transform unit

圖4 單元變換2Fig.4 2 transform unit

故Delaunay三角單元有限元求解橢圓型偏微分方程的流程圖如圖5.

圖5 有限元方法求解流程圖Fig.5 Finite element method for solving the flow chart

3 數(shù)值例子

本文有限元方法采用的是Lagrange二次元插值型，在計(jì)算數(shù)值積分時(shí)可用到了高斯積分技術(shù)[1]，所有的程序是在Matlab 2012編程實(shí)現(xiàn)的.本文單元?jiǎng)澐植介L取的是0.1.通過Matlab程序運(yùn)行，得到真實(shí)值和半有限元法解如圖6.

圖6 有限元求解對(duì)比圖Fig.6 Finite element solution comparison chart

上圖5中，紅色線表示真實(shí)解，綠色線表示有限元數(shù)值解.左邊圖，在求有限元數(shù)值解時(shí)，沒有進(jìn)行Delaunay三角劃分和單元自適應(yīng)編號(hào)處理.右邊圖，在求有限元數(shù)值解時(shí)，單元?jiǎng)澐謺r(shí)，對(duì)三角單元進(jìn)行了Delaunay三角化處理，并且對(duì)單元節(jié)點(diǎn)采用了自適應(yīng)編號(hào)排序處理.從上面兩圖比較可以看出，通過對(duì)三角單元進(jìn)行了Delaunay三角化處理和自適應(yīng)編號(hào)處理，可以減少舍入誤差，計(jì)算結(jié)果的精度更高.從右圖可以得出，紅色線和綠色線十分吻合，可以說明本文有限元方法得到的數(shù)值解，逼近真實(shí)解，即本文的有限元方法是可行的.

4 結(jié)論

有限元方法作用于微分方程，其本質(zhì)是通過數(shù)學(xué)方法，將微分方程轉(zhuǎn)化成積分方程，對(duì)問題區(qū)域進(jìn)行單元?jiǎng)澐郑卧逯岛瘮?shù)，通過在各個(gè)單元上運(yùn)用有限元方法，應(yīng)用高斯數(shù)值積分求取單元上的系數(shù)矩陣，通過單元關(guān)系矩陣將單元上的系數(shù)矩陣，組裝成總系數(shù)矩陣，求解有限元方程，得到偏微分方程的有限元方法的數(shù)值解.本文將Delaunay三角化運(yùn)用到三角單元?jiǎng)澐郑?duì)單元節(jié)點(diǎn)采用自適應(yīng)編號(hào)排序，減少了舍入誤差，從而提高了數(shù)值解的精度.

[1] 倪 興. 常微分方程數(shù)值解法及其應(yīng)用[D]. 合肥：中國科學(xué)技術(shù)大學(xué), 2010．

[2] GIRAULT V, RAVIART P A. Finite element method for Navier-Stokes equations: theory and algorithms[M]. Berlin: Heidelber Springer-Verlag, 1981．

[3] BLACK F, SCHOLES M. The pricing of options and corporate liabilities[ J] . Journal of Political Economy , 1973 , 81(4): 633-654．

[4] Ho-Le K. Finite element mesh generation methods: a review and classification[J]. Computer Aided Design, 1988, 20(1):27-38．

[5] 傅永華.有限元分析基礎(chǔ)[M]. 武漢:武漢大學(xué)出版社, 2003．

[6] Zienkeiwicz O C. The finite element method[M]. 3rd edition, McGraw-Hill, 1977．

[7] COGGON J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(2):132-155．

[8] MILLER G L, TALMOR D, TENG S H, et al. Control valume meshes using sphere packing: generation, refinement and coarsening[J]. Proceeding of 5th International Meshing Roundtable, 1996, 47-61．

[9] Green P J, Sibson R. Computing Dirichlet tessellation in the plane[J]. The Computer Journal, 1987, 2(2):168-173．

[10] Bowyer A. Computing Dirichlet Tessellations[J]. The Computer Journal, 1981, 24(2):162-166．

[11] 楊 欽. 限定Delaunay三角網(wǎng)絡(luò)剖分技術(shù)[M].北京: 電子工業(yè)出版社， 2005.

[12] 林 群. 微分方程數(shù)值解法基礎(chǔ)教程[M].第二版.北京:科學(xué)出版社， 2003．

[13] 華東師范大學(xué)數(shù)學(xué)系. 數(shù)學(xué)分析[M]. 北京: 高等教育出版社, 2001．

[14] 羅智中. 一種任意二維區(qū)域圖形的三角剖分方法[J].機(jī)床與液壓, 2002(6):145-146．

Triangle finite element numerical solution of elliptic partial differential equations

YU Hui1, YANG Zhenhao1, FENG Jiming1, QU Shaocheng2

(1.College of Science，China Three Gorges University，Yichang, Hubei 443000;2.Department of Electronics and Information Engineering，Central China Normal University，Wuhan 430079)

In the present work, the finite element method for triangular element is studied based on the elliptic partial differential equation. With the development of computer graphics, the triangle element division has made great achievements. High-quality triangular elements are able to be generated, which improved the accuracy of numerical solutions for partial differential equations. Here, the triangular element satisfying the Delaunay conditions is calculated through unit division of the region with the node incremental algorithm. All nodes of the triangular element are adopted with adaptive numbers and the numerical solution of elliptic partial differential equations is obtained utilizing the finite element method of triangular elements. Finally, through numerical experiments, the finite element method is compared with the traditional one. Results indicate that the finite element method of triangular element reduces the rounding error and elevates the calculation precision.

delaunay triangular element; adaptive number; rounding error

2016-03-16.

國家自然科學(xué)基金項(xiàng)目(61273183, 61374028, 61304162).

1000-1190(2016)04-0489-07

O175.25

A

*E-mail: yuhui@ctgu.edu.cn.