微分方程有限元法的邊值與初值問題相容性處理

李純清， 龔興厚， 張高文

(湖北工業(yè)大學材料科學與工程學院， 湖北 武漢 430068)

有限元法是一種離散的變分法或基于變分原理的離散方法．它吸收了古典Ritz-Galerkin方法和有限差分法的精髓，既以變分原理為基礎，又對解空間進行離散近似，同時還利用分段插值多項式特點，從而克服了微分方程數(shù)值求解許多刺手問題1-2．

1 數(shù)值求解中的邊值及初值問題處理原理與方法

以比較典型的一類微分方程及其各類邊界條件為例[3]:

)+q(x)y=

f(x),x∈(a,b);

(1)

p(x)∈C1[a,b],q(x),f(x)∈C[a,b].

式(1)的Galerkin變分形式略推

(Ly,v)=(f,v),

即

(2)

對式(2)左邊第一項進行分部積分：

v(a)ga-v(a)αay(a)-v(b)gb+v(b)αby(b)+

將分部積分結(jié)果代入式(2)并整理后得

)dx-

v(a)αay(a)+v(b)αby(b)=

(3)

將上面等式中兩個積分項中y,v進行有限維化，也就是構(gòu)造連續(xù)函數(shù)集合S的有限N+1維子集SN+1，使得y,v∈SN+1=span{φ0,φ1,φ2,…,φn}，記為yN+1,vN+1(稱為Galerkin連續(xù)變分問題離散近似化)，并將積分區(qū)間[a,b]分段化a,x1,x2,…,xn-1,b，在每個結(jié)點上定義插值基函數(shù)φ0,φ1,…,φn(即有限元化)，得到所謂總體剛度矩陣和總體載荷向量間的方程，取插值基函數(shù)φi為一次項式，其一般形式如

結(jié)合式(3)，可將式(1)中的邊值及初值條件的處理

(4)

2 實例驗證

例 1： 初值問題

其準確解為y(x)=ex-e-x+1.

將求解區(qū)間分段化，取6個節(jié)點分別為0,0.2,0.4,0.6,0.8,1，離散為5個子區(qū)間[0,0.2]、[0.2,0.4]，[0.4,0.6]，[0.6,0.8]，[0.8,1]，在各個區(qū)間的節(jié)點上取一次插值基函數(shù)，用MatLab容易實現(xiàn)例1對應的Galerkin變分形式的線性方程組為

對例1給定初值條件進行處理

1×2+α0×1=g0, 1×dy5+α1×y5=g1

令α0=0;則g0=2

令α1=0;則g1=dy5

將α0=0,g0=2 ,α1=0,g1=dy5(待定)按式(4)處理方法可得到例1的初值處理后的線性方程組

對上述方程組利用MatLab的符號運算功能易求取其解

其中已知y0=y(0)=1，即-1.6193+0.8476×dy5=1，可得到dy5=3.0901，并代入上式得到例1在節(jié)點0,0.2,0.4,0.6,0.8,1的數(shù)值解(表1)，與理論解十分接近．

表1 例1數(shù)值解

本例也可令α0=2,則g0=4; 令α1=1,則g1=dy5+y5．但需要建立方程組求出待定的dy5與y5，具體過程

(5)

建立方程組求解dy5,y5

解之得

代入式(5)，也可求出在節(jié)點0,0.2,0.4,0.6,0.8,1的數(shù)值解，結(jié)果與前述方法基本相同．

例2：兩點邊值問題：

y'(1)=1,y'(2)=4

用MatLab可求取其理論解表達式，取x={1，1.25，1.5，1.75，2}處的理論值，則y(x)分別為7.6515，8.0707，8.7289，9.5469，10.4904.

同理，插值基函數(shù)取二次項式，用MatLab容易實現(xiàn)例2對應的Galerkin變分形式的線性方程組為

對例2給定邊值條件進行處理如

1×1+α1×y0=g1

4×4+α2×y4=g2

令α1=0;則g1=1

令α2=0;則g2=16

將α1=0,g1=1 ,α2=0,g2=16按式(4)處理方法可得到例2的邊值問題處理后的線性方程組

其解與理論解比較非常接近(表2)．

表2 例2數(shù)值解

例3：兩點邊值問題：

y=-xy(1)=1,y(2)=4.

其理論解在x為1，1.25，1.5，1.75，2值時，y(x)依次對應為1.0000，1.9435，2.6945，3.3626，4.0000.

對例3給定邊值條件進行處理

1×dy0+α1×1=g1

4×dy2+α2×4=g2

令α1=1,則g1=1+dy0

令α2=2,則g2=8+4×dy4.

將α1=1,g1=1+dy0;α2=2,g2=8+4×dy4(其中dy0,dy4待求)，按式(4)處理方法可得到例3的邊值問題處理后的線性方程組

同樣，借助MatLab符號運算功能，求出上述方程的解為

(6)

由于y0=y(1)=1,y4=y(2)=4,可如下建立方程組，求出dy0,dy4．

再將dy0,dy4結(jié)果代入式(6)，可計算出例3在節(jié)點x=1，1.25，1.5，1.75，2處的數(shù)值解(表3)，與理論值也非常接近．

表3 例3數(shù)值解

3 討論與總結(jié)

1)無論是微分方程的初始條件，還是邊界條件，本處理方法都適用．尤其是非齊次邊界條件，比打靶法優(yōu)越．

2)本文提出的對微分方程的初值和邊值問題的處理方法與技巧主要是針對它的Galerkin變分形式進行討論的，但也同樣適用于微分方程的Ritz變分形式，因為Ritz變分形式只是Galerkin變分形式的一種特殊情況．

3)處理的依據(jù)是微分方程對應的等價泛函變分形式．在尋求變分形式過程中，微分方程特定的初始條件和(或)邊界條件的也映射為泛函變分形式另一種特定形式，第三類邊界條件可以說是兩者間的映射公式．

4)一般而言，并不是所有的邊值與初值問題存在對應的變分形式，這就使得有限元法應用受到限制．其實，如果給定的邊值或初值條件能與微分方程相兼容，認為有限元法就可應用，比如實際物理問題(溫度場，應變場等)的微分方程描述．

5)將微分方程不進行變分方法，而直接應用有限元法進行數(shù)值求解，此問題有待探索．因為有限元思想是獨立，當與Ritz-Galerkin變分近似解法相結(jié)合，才有用有限元法數(shù)值求解微分方程．

[參考文獻]

[1] 施妙根. 科學和工程計算基礎[M]. 北京：清華大學出版社, 2002(7)：239-266.

[2] 老大中. 變分法基礎[M]. 北京：國防工業(yè)出版社, 2004(9):285-323.

[3] 姚端正. 數(shù)學物理方法[M]. 北京：科學出版社, 2010(3):(220-235.