蔡姝婷

（福建江夏學(xué)院 數(shù)理教研部，福建 福州 350108）

偏微分方程解的數(shù)值驗(yàn)算

蔡姝婷

（福建江夏學(xué)院 數(shù)理教研部，福建 福州 350108）

本文提出一種偏微分方程解的數(shù)值驗(yàn)算方法.首先應(yīng)用有限元的方法求偏微分方程的數(shù)值解，然后將偏微分方程的解轉(zhuǎn)化為一個緊算子的不動點(diǎn)，我們驗(yàn)證算子在某一個集合中滿足不動點(diǎn)定理，從而偏微分方程的解存在.接著通過分析的方法將不動點(diǎn)定理的條件轉(zhuǎn)化為計算機(jī)可以計算的條件，最后結(jié)合區(qū)間驗(yàn)算的方法避免計算機(jī)計算過程中存在的截斷誤差，從而證明所求近似解的附近存在精確解.

偏微分方程；數(shù)值驗(yàn)算；不動點(diǎn)定理；區(qū)間驗(yàn)算

1 引言

偏微分方程在自然科學(xué)和技術(shù)科學(xué)領(lǐng)域中都具有很廣泛的應(yīng)用.在每年的各類型數(shù)學(xué)建模競賽中，運(yùn)用偏微分方程進(jìn)行建模是一種常見的建模方法.如反應(yīng)擴(kuò)散型方程中的食餌－捕食者模型，可以模擬兩個不同的種族間的相互作用，其中一個種族捕食另一個種族；Schnakenberg模型描述發(fā)育生物學(xué)中的一些模式形成等等.如何求這些偏微分方程的解是能否完成一個數(shù)學(xué)建模過程的重要的問題.到目前為止，能夠運(yùn)用解析的方法求解的偏微分方程類型相對有限.但是，隨著計算機(jī)計算能力的增強(qiáng)以及數(shù)值計算方法的發(fā)展，偏微分方程的數(shù)值解法越來越多，通常運(yùn)用一些有限元的方法將偏微分方程離散化后進(jìn)行求解，之后需要求出相應(yīng)的誤差.有部分的數(shù)學(xué)軟件如Matlab可以用自帶的軟件包直接進(jìn)行數(shù)值計算，但這樣的數(shù)值計算過程都帶有一定的誤差，無法嚴(yán)格地證明方程的解確實(shí)在某個范圍內(nèi)是存在的.本文的方法是對傳統(tǒng)的數(shù)值計算方法的發(fā)展，提出通過數(shù)值驗(yàn)算的方法驗(yàn)算偏微分方程解的存在性.我們先用有限元的方法求出偏微分方程的數(shù)值解，然后將偏微分方程的解轉(zhuǎn)化為一個方程的不動點(diǎn)，根據(jù)不動點(diǎn)定理，驗(yàn)證在一個小區(qū)域內(nèi)，這個方程的不動點(diǎn)存在.這個驗(yàn)證的過程是通過計算機(jī)來完成的，而計算機(jī)在內(nèi)部計算的過程中，存在一定的截斷誤差，我們應(yīng)用區(qū)間演算[1]的方法來避免這個誤差.整個計算過程我們將在Matlab上實(shí)現(xiàn).通過我們的方法，一方面可以得出方程解的大致形態(tài)，另一方面，也可以用計算機(jī)進(jìn)行嚴(yán)格的證明解確實(shí)是存在的以及解所存在的小范圍.

我們選取較為簡單的偏微分方程說明本文的方法.考慮具有Dirichlet邊界條件的偏微分方程

其中Ω=(0,1)×(0,1),f(x,y)是關(guān)于x,y的函數(shù)，d是一個常數(shù).

方程(1)近似解的求法有很多種，在[2]中介紹了多種求解偏微分方程數(shù)值解的有限元方法.本文同樣將應(yīng)用有限元的方法求解方程(1),只是對于基函數(shù)的選取有所不同.

首先，對于任意的整數(shù)m，我們定義Hm=Hm(Ω)為m階的L2-Sobolev空間.同時我們定義H01:=H01(Ω):={u∈H1|u=0 on Ω}，并具有內(nèi)積〈?,ψ〉=(▽?,▽ψ)，其中(·,·)表示L2(Ω)空間上的內(nèi)積.我們考慮方程（1）在H01空間中解的存在性.

我們選取基函數(shù)φi1i2為

對于一個固定的非負(fù)整數(shù)N，我們按以下規(guī)則對基函數(shù)重新排序，

從而基函數(shù)成為

假設(shè)H01中的所有元素ψ可展開為傅利葉級數(shù)

選取適當(dāng)?shù)腘，我們在有限維空間

2 近似解

在這一節(jié)中，我們將在空間SN中計算方程(1)的近似解.對任意的vN∈SN令

因此我們現(xiàn)在的問題是，求uh∈Sh,使得

將上述問題寫成矩陣形式，即

其中

我們運(yùn)用Matlab即可求出方程(1)的數(shù)值解.

3 解的數(shù)值驗(yàn)算

我們知道對于任意的ψ∈L2(Ω)下列的泊松方程存在一個唯一的根?∈H1∩H'0：

將方程(2)的解記為?=L-1ψ則L-1:L2→H01是一個緊算子.

令K(u):=L-1(f(x,y)+du),那么K是H01上的緊算子，并且方程(1)的解為算子K的不動點(diǎn)：

因此由Schauder不動點(diǎn)定理，如果我們能在H01中找到一個非空的，有界的，凸的閉集合滿足

那么在K(U)中就存在一個元素u使得u=Ku.我們的做法是在計算機(jī)中構(gòu)造一個侯補(bǔ)的集合U=UN⊕U⊥,其中UN?SN,U⊥?SN⊥.這里SN⊥表示SN在H01空間中的正交補(bǔ)空間.這樣，我們可以將驗(yàn)證條件(4)轉(zhuǎn)化為

Moore在[1]中介紹了區(qū)間演算的方法.由于在計算機(jī)計算的過程中存在截斷誤差，比如在數(shù)學(xué)軟件Matlab中，以雙精度符點(diǎn)型數(shù)據(jù)進(jìn)行計算，當(dāng)我們計算1/3時，計算機(jī)內(nèi)部以64位小數(shù)進(jìn)行計算，截斷了小數(shù)點(diǎn)后的值，由此產(chǎn)生了誤差.我們應(yīng)用區(qū)間演算來避免這種誤差，如同樣計算1/3時，計算機(jī)內(nèi)部以[0.33333333333333,0.33333333333334]計算.此時既有1/3落于這個小區(qū)間內(nèi)，避免了截斷誤差，同時這個區(qū)間的范圍較小，滿足一定的精度要求.為了計算的順利進(jìn)行，我們定義了區(qū)間演算的加減乘除法.對于實(shí)數(shù)集上的任意兩個區(qū)間，X=[a,b],Y=[c,d],我們定義下列運(yùn)算：

其中x=min{ac,ad,bc,bd},y=max{ac,ad,bc,bd}.同時還定義

例如：[1,2]·[-2,1]+[0,1]=[-4,3].另外，區(qū)間X的幅度定義為：d(X)≡b-a.絕對值為：|X|≡max{|a|,|b|}以及中點(diǎn)為m(X)≡(a+b)/2.Rump根據(jù)這些運(yùn)算規(guī)則，在[3]中建立了一個基于區(qū)間演算的區(qū)間實(shí)驗(yàn)室，實(shí)現(xiàn)了在計算機(jī)上完成區(qū)間演算的過程.

我們這里將用這種區(qū)間演算的方法驗(yàn)證條件(4)的成立.對于條件(4)的有限維部分PNK(U)?UN的驗(yàn)證，我們把候補(bǔ)集合UN定義為

其中Bi是由方程

所確定的.

而對于無限維部分的驗(yàn)證，我們需要一個引理：

引理1 對任何的?∈H2∩H01,存在一個正常數(shù)CN，使得

定理1如果

成立，那么在U中存在一個方程K(u)=u的解.

根據(jù)這個定理，我們可以得到驗(yàn)證解存在性的一個算法：

算法1

4 數(shù)值結(jié)果

我們選取具體的d及f進(jìn)行計算.若d=1,f(x,y)=x+y.取N=20.我們用Matlab軟件得到方程(1)的近似解如圖1所示.接著利用區(qū)間演算實(shí)驗(yàn)室[3]以及算法1進(jìn)行數(shù)值驗(yàn)算，得到區(qū)間長度的最大值為0.004,β0=0.0175.本文所得到的計算結(jié)果是在HP ENVY 6 Notebook PC [Intel(R)Core(TM)i5-3317U@1.70 GHz]電腦上用Matlab (Ver.7.0)完成的.

圖1

注：本文的方法針對于相對較為簡單的方程(1)進(jìn)行求解，具有一定的局限性.如果所遇到的方程關(guān)于u是非線性的，同樣是可以用數(shù)值驗(yàn)算的方法進(jìn)行求解.在求近似解時，可以先用Newton法，而后的數(shù)值驗(yàn)算過程可以采用一種近似的牛頓法，結(jié)合不動點(diǎn)定理加以驗(yàn)證.這是我們今后的研究方向.

〔1〕Moore R.E.Intervalanalysis[M].Englewood Cliffs: Prentice-Hall,1966.

〔2〕馬昌鳳.有限元方法講義[Z].2006.

〔2〕Rump S.M.,INTLAB-INTerval LABoratory[Z],a Matlab toolbox for verified computations,Version 5.5.Inst,Informatic,TU Hamburg-Harburg,http://www.ti3.tuharburg.de/rump/intlab/.

O175.2

A

1673-260X（2015）09-0007-03

本文由教育部留學(xué)回國人員科研啟動基金資助項目資助