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

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

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

?

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

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

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

本文以橢圓型偏微分方程問題為背景研究三角形單元的有限元方法.隨著計算機圖形學(xué)的發(fā)展，三角形單元劃分取得了巨大的成就，可以得到質(zhì)量非常好的三角形單元，進而提高偏微分方程的有限元方法的數(shù)值解的精度.本文采用節(jié)點增量算法，對問題區(qū)域進行三角形單元劃分，得到的三角形單元滿足Delaunay條件，再對三角形單元的所有節(jié)點采用自適應(yīng)編號，最后運用三角形單元的有限元方法得到橢圓型偏微分方程的數(shù)值解.通過數(shù)值實驗，得出相比傳統(tǒng)的三角形單元的有限元方法，本文的三角形單元的有限元方法減小了舍入誤差，提高了計算精度.

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

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

有限元方法[4]是求解偏微分方程數(shù)值解的重要方法.有限元求解偏微分方程時，一個重要的步驟是[5]，運用數(shù)學(xué)知識將偏微分方程轉(zhuǎn)化為積分方程.對偏微分方程的原問題區(qū)域進行單元劃分[6-8]，如四邊形單元劃分和三角單元劃分.20世紀80年代三角網(wǎng)格剖分在眾多領(lǐng)域引起了廣泛的關(guān)注.英國Bath大學(xué)數(shù)學(xué)分校的P.J Green和R. Sibson等從數(shù)學(xué)的角度對三角網(wǎng)格剖分進行了研究[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)格.

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

1 有限元構(gòu)思

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

(1)

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

.

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

(2)

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

(3)

其中

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

(4)

記

于是得到如下式子

(5)

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

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

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

(6)

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

從問題的積分形式可以看出，二階橢圓型偏微分方程的條件是u的二階偏導(dǎo)存在，通過運用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的多項式空間，則本文所要用到的多項式樣條函數(shù)空間，可以寫為

(7)

(8)

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

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

2.1三角形單元劃分

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

圖1 圓形區(qū)域三角形單元劃分Fig.1 Circular area triangle unit division

2.2單元變換

圖2 單元變換Fig.2 Unit transform

設(shè)

2.3單元插值函數(shù)

設(shè)

運用Matlab編程計算，可得：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.

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

Npup+Nquq+Nrur.

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

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

因為Lagrange二次型插值函數(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)編號排序

下圖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二次元插值型，在計算數(shù)值積分時可用到了高斯積分技術(shù)[1]，所有的程序是在Matlab 2012編程實現(xiàn)的.本文單元劃分步長取的是0.1.通過Matlab程序運行，得到真實值和半有限元法解如圖6.

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

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

4 結(jié)論

有限元方法作用于微分方程，其本質(zhì)是通過數(shù)學(xué)方法，將微分方程轉(zhuǎn)化成積分方程，對問題區(qū)域進行單元劃分，建立單元插值函數(shù)，通過在各個單元上運用有限元方法，應(yīng)用高斯數(shù)值積分求取單元上的系數(shù)矩陣，通過單元關(guān)系矩陣將單元上的系數(shù)矩陣，組裝成總系數(shù)矩陣，求解有限元方程，得到偏微分方程的有限元方法的數(shù)值解.本文將Delaunay三角化運用到三角單元劃分，并對單元節(jié)點采用自適應(yīng)編號排序，減少了舍入誤差，從而提高了數(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].機床與液壓, 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é)基金項目(61273183, 61374028, 61304162).

1000-1190(2016)04-0489-07

O175.25

A

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