兩點(diǎn)邊值問題的三類邊值條件的有限元解法實(shí)現(xiàn)

盧仁洋 于陸洋 江山

摘 要：文章研究二階微分方程的兩點(diǎn)邊值問題，使用有限元方法對(duì)三類不同的邊值條件具體進(jìn)行討論和處理。對(duì)于可齊次化的Dirichlet、Neumann邊值， 給出相應(yīng)分析以簡(jiǎn)化和規(guī)范計(jì)算步驟。對(duì)于Robin邊值， 基于之前的分析給出實(shí)現(xiàn)技巧以達(dá)到有效的數(shù)值模擬。

關(guān)鍵詞：有限元方法；Dirichlet邊值；Neumann邊值；Robin邊值

中圖分類號(hào)：O172 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：2096-000X（2018）05-0058-03

Abstract： A two-point boundary value problem of the second-order ordinary differential equation is studied in this paper， and a finite element method is used to deal with three kinds of boundary conditions in detail. The corresponding analysis is provided to simplify and standardize the steps for the homogeneous Dirichlet， Neumann boundary values. For the Robin boundary value， we utilize the previous analyses to present a specific strategy for a good performance in numerical simulations.

Keywords： finite element method； Dirichlet boundary； Neumann boundary； Robin boundary

一、本文模型及介紹

二階微分方程的兩點(diǎn)邊值問題是科學(xué)工程計(jì)算中的經(jīng)典問題，也是微分?jǐn)?shù)值解法的必要基礎(chǔ)[1，2]。常見的做法為有限差分法[3，4]、有限元方法[5-7]、有限體積法[8，9]等。

本文考慮兩點(diǎn)邊值問題（BVP）

其中變系數(shù) 。

上述區(qū)間I也可以是I=[a，b]，但為研究方便，能通過變換x'=轉(zhuǎn)化為I=[0，1]。

針對(duì)上述方程，根據(jù)A1，A2，B1，B2取值的不同，分成三類邊值條件：

1. Dirichlet邊值：A1≠0，B1≠0，A2=0，B2=0，即u（0）=？琢/A1，u（1）=？茁/B1。

2. Neumann邊值：A1=0，B1=0，A2≠0，B2≠0，即u'（0）=？琢/A2，u（1）=？茁/B2。

3. Robin邊值：A1≠0，B1≠0，A2≠0，B2≠0。

本文采用有限元法（Finite Element Method， FEM）處理模型，首先需確定變分形式。對(duì)（1.1）移項(xiàng)后Lu-f=0，兩邊同乘v并在[0，1]上積分，得：

即有： （1.2）

其中a（u，v）=（p+qv+ruv）dx，（f，v）=fvdx，H1為一階可導(dǎo)的完全內(nèi)積空間。設(shè)Uh是U的n維子空間，選取Un的一組基底？漬1，？漬2…，？漬n作為基函數(shù)，有限元解即為uh=u？漬，取線性元基函數(shù)？漬1，？漬2…，？漬n計(jì)算。通過||u-uh||∞計(jì)算誤差、log2（）計(jì)算收斂率，其中u真解、uh為有限元解。

二、三類邊值問題的分析和算例

根據(jù)邊值條件是否齊次，又可具體細(xì)分。

（一）Dirichlet邊值

（1）？琢=？茁=0，即u（0）=0，u（1）=0

（2）？琢=0，？茁≠0，即u（0）=0，u（1）=？茁/B1

處理（2）的做法：令w（x）=u（x）-u（1），則有：

即 ，就轉(zhuǎn)化為（1）的情形。

（3）？琢≠0，？茁=0，即u（0）≠？琢/A1=？琢D，u（1）=0?？梢酝ㄟ^令：（x）=u（x）+u（0），后續(xù)處理與（2）一致。

（4）？琢≠0，？茁≠0，即u（0）≠？琢/A1=？琢D，u（1）=？茁/B1=？茁D，可齊次化左邊界，令w（x）=u（x）-u（0），則有：

即 ，進(jìn)而轉(zhuǎn)化為（2）。

算例1：式（1.1）中，取p（x）=x q（x）=3x，r（x）=-2x，f（x）=-2e2x-ex。

A1=1，A2=0，？琢=1，B1=1，B2=0，？茁=e2+e精確解為u=e2x+xex。

通過上述分析，令w（x）=u（x）-u（0）-（u（1）-u（0））=u（x）-（e2+e-1）x-1， 則有

經(jīng)典有限元法求出wh后，得到uh=wh+（e2+e-1）x+1。 計(jì)算結(jié)果如下圖表：

可以看出，在N=32較粗網(wǎng)格即有10-3級(jí)誤差，隨著網(wǎng)格不斷剖分、有限元法的逼近效果越來越好，并且收斂率達(dá)到穩(wěn)定的2階收斂，這是符合預(yù)期結(jié)果的。

（二）Neumann邊值

（1）？琢=？茁=0，即u'（0）=0，u'（1）=0。

（2）？琢=0，？茁≠0，即u'（0）=0，u'（1）=？茁/B2=？茁N。

（3）？琢≠0，？茁=0，即u'（0）≠？琢/A2=？琢N，u'（1）=0。

（4）？琢≠0，？茁≠0，即u'（0）≠？琢/A2=？琢N，u'（1）=？茁/B2=？茁N。

通過與Dirichlet邊值類似的處理方法，先齊次化左邊界，再齊次化右邊界。具體做法：令w（x）=u（x）-u'（0）x-■ （u'（1）-u'（0））x2

即有

算例2：式（1.1）中，取p（x）=x q（x）=3x，r（x）=-2x，f（x）=-2e2x-ex。

A1=0，A2=1，？琢=3，B1=0，B2=1，？茁=2e2+2e該方程的精確解為u=e2x+xex。

由上述分析，令w（x）=u（x）-3x-■（2e2+2e-3）x2，原問題轉(zhuǎn)化為

經(jīng)典有限元法求出vh后，得到uh=wh+■（2e2+2e-3）x2+3x。

表2 Neumann邊值問題的誤差與收斂率

（三）Robin邊值

u'（0）=（？琢-A1u（0））/A2=？琢R1-？琢R2u（0），u'（1）=（？茁-B1u（1）））/B2=？茁R1-？茁R2u（1）

由于Robin邊值特性，導(dǎo)數(shù)項(xiàng)與函數(shù)項(xiàng)混合，無法使用齊次化方法。這里給出一般化的解決方法：在形成總剛度矩陣以及右端載荷向量時(shí)，根據(jù)邊值條件相對(duì)應(yīng)地修改矩陣和向量，從線性元性質(zhì)出發(fā)，實(shí)際需要修改的是矩陣的第1×1項(xiàng)和第（N+1）×（N+1）項(xiàng)以及載荷向量的第1項(xiàng)和第N+1項(xiàng)。

算例3：式（1.1）中，取p（x）=x，b（x）=3x，c（x）=-2x，f（x）=-2e2x-ex。

A1=1，A2=1，？琢=4，B1=-2，B2=1，？茁=0該方程的精確解為 u=e2x+xex。

利用經(jīng)典有限元法形成總剛度矩陣A以及載荷向量 F，構(gòu)成了（N+1）×（N+1）階數(shù)，將Robin邊值的導(dǎo)數(shù)項(xiàng)用函數(shù)項(xiàng)表示并代入（1.2）式得：

表3 Robin邊值問題的誤差與收斂率

通過詳細(xì)分析和具體驗(yàn)證，我們驗(yàn)證了三類邊值中無論哪類模型，使用有限元法和數(shù)學(xué)技巧處理都得到了有效精度和穩(wěn)定收斂的2階結(jié)果。本文是作者學(xué)習(xí)過程中對(duì)于有限元知識(shí)的一些總結(jié)和融合，將三類邊值問題做了具體的分析以及數(shù)值案例的展示，也對(duì)線性有限元的誤差進(jìn)行了實(shí)際驗(yàn)證。

參考文獻(xiàn)：

[1]李榮華.偏微分方程數(shù)值解法[M].北京：高等教育出版社，2012.

[2]陸金甫.偏微分方程數(shù)值解法[M].北京：清華大學(xué)出版社，2016.

[3]許芝卉，王鳳艷.用差分方法求解兩點(diǎn)邊值問題[J].山西大同大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，25（4）：10-11.

[4]趙秉新，鄭來運(yùn).兩點(diǎn)邊值問題的差分方法及快速求解[J].華南師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2010（4）：27-30.

[5]Huang Z Y. Tailored finite point method for the interface problem[J].Networks and Heterogeneous Media，2009，4（1）：91-106.

[6]付小龍，李名，李小瑞.兩點(diǎn)邊值問題的有限元算法[C].北京力學(xué)會(huì)學(xué)術(shù)年會(huì)，2014.

[7]Constantinou P， Xenophontos C. Finite element analysis of an exponentially graded mesh for singularly perturbed problems[J]. Computational Methods in Applied Mathematics，2015，15（2）： 135-143.

[8]Cao W， Zhang Z， Zou Q. Superconvergence of any order finite volume schemes for 1D general elliptic equations[J].Journal of Scientific Computing，2013，56（3）：566-590.

[9]魏保軍，張武軍，石金娥.兩點(diǎn)邊值問題有限體積法的一種加權(quán)模估計(jì)[J].廣西師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2015，33（3）：75-78.