張靜, 王彩華

(天津師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院,天津 300387)

1.引言

奇異擾動(dòng)方程一般是指微分方程中最高階導(dǎo)數(shù)項(xiàng)含一個(gè)小參數(shù), 此類問(wèn)題廣泛出現(xiàn)在應(yīng)用數(shù)學(xué)和工程學(xué)的許多領(lǐng)域.求解奇異擾動(dòng)方程的方法有很多, 如有限差分法[1-2]、有限元法[3-4]、有限體積法[5-6]等.本文研究一種奇異擾動(dòng)兩點(diǎn)邊值問(wèn)題: 含參數(shù)ε的含源反應(yīng)擴(kuò)散問(wèn)題, 其微分方程模型為

針對(duì)一維對(duì)流擴(kuò)散反應(yīng)方程的數(shù)值解有較多研究, 如文[7]根據(jù)一階和二階導(dǎo)數(shù)的padé型緊致差分格式推導(dǎo)出一種四階格式; 文[8]首先研究常系數(shù)情形下的差分格式, 然后在此基:上修正, 提出了變系數(shù)對(duì)流擴(kuò)散反應(yīng)方程的指數(shù)型高精度差分格式; 文[9]引入了指數(shù)擬合因子,結(jié)合等距網(wǎng)格給出了一種有限差分格式, 適合于求解具有雙邊界層的兩點(diǎn)邊值問(wèn)題.

針對(duì)反應(yīng)擴(kuò)散模型(1.1), 文[10]設(shè)計(jì)了一種六階三點(diǎn)差分格式FD1[10]:

其中

文[11]針對(duì)對(duì)流擴(kuò)散反應(yīng)方程提出了差分格式FD2[11], 應(yīng)用于模型方程(1.1)得到如(1.2)的三點(diǎn)格式, 其中

該格式具有四階精度.本文實(shí)驗(yàn)部分主要與文[10-11]給出的差分格式進(jìn)行了數(shù)值比對(duì).

雖然奇異擾動(dòng)問(wèn)題已有較多數(shù)值方法, 但一些方法適用于參數(shù)較大的情形, 當(dāng)邊界層非常窄時(shí)失效, 能否設(shè)計(jì)出精度更高, 適用性更強(qiáng)的數(shù)值格式仍值得進(jìn)一步研究.文[12]針對(duì)無(wú)源對(duì)流擴(kuò)散問(wèn)題設(shè)計(jì)了兩類修正差分格式: 其橫向系列修正格式結(jié)合等距網(wǎng)格能高精度求解參數(shù)不是很小的對(duì)流擴(kuò)散問(wèn)題, 精度可分別達(dá)到二階、四階、六階和八階, 其縱向格式能高精度求解參數(shù)很小時(shí)的對(duì)流擴(kuò)散問(wèn)題, 數(shù)值實(shí)驗(yàn)表明, 文[12]提出的兩類方法與八種參考格式比對(duì)數(shù)值結(jié)果更好.本文擬沿用文[12]的差分格式設(shè)計(jì)思路, 展開(kāi)對(duì)反應(yīng)擴(kuò)散問(wèn)題數(shù)值解的研究.

為表述簡(jiǎn)單, 下文的系列橫向差分格式簡(jiǎn)記為HDS1, HDS2, HDS3 (HDS: Horizontal Difference Scheme)等; 縱向差分格式記為VDS1, VDS2 (VDS: Vertical Difference Scheme)等.

2.橫向與縱向系列修正差分格式的設(shè)計(jì)思路

設(shè)函數(shù)u(x)充分光滑, 當(dāng)i= 1,2,··· ,N -1時(shí), 三點(diǎn)模板上將ui-1、ui+1在點(diǎn)xi處進(jìn)行Taylor展開(kāi)

于是, 當(dāng)n=2,3,···時(shí), 遞推公式為

依據(jù)遞推公式(2.5), 可列出u(x)的部分高階導(dǎo)數(shù)降階組式如下

橫向系列差分格式設(shè)計(jì)思路: 當(dāng)ε不是很小時(shí),把(2.6)中部分降階式代入泰勒展式(2.1),(2.2)中, 關(guān)于h的冪次保留的項(xiàng)數(shù)不同則截?cái)嗾`差不同.例如, 式(2.1),(2.2)只取到u′′(x)項(xiàng), 整理后可得到二階差分格式HDS1(即中心差分格式CDS); 式(2.1),(2.2)如保留到u(3)(x),u(4)(x)項(xiàng),整理后可得四階差分格式HDS2; 繼續(xù)修正可得六階差分格式HDS3.

縱向系列差分格式設(shè)計(jì)思路: 當(dāng)ε非常小時(shí), 縱向觀察組式(2.6), 考慮含1/ε的負(fù)指數(shù)冪次越大的項(xiàng)其絕對(duì)值可能越大.如組式(2.6)僅保留縱向關(guān)于1/ε指數(shù)冪次最大項(xiàng)代入泰勒展式(2.1),(2.2), 舍去其它項(xiàng), 整理可得縱向差分格式VDS1; 在VDS1的基:上繼續(xù)修正,保留1/ε指數(shù)冪次最大項(xiàng)及次大項(xiàng), 整理可得縱向差分格式VDS2; 類似地, 進(jìn)一步修正可得VDS3等.

如需進(jìn)一步提高數(shù)值解的精度, 且不考慮計(jì)算機(jī)舍入誤差, 橫向、縱向差分格式可一直修正下去.

3.橫向系列修正差分格式

Ⅰ 二階橫向差分格式HDS1(中心差分格式CDS)

由帶微分余項(xiàng)的Taylor展式知

其中ξi ∈(xi-1,xi+1).將式(3.3)代入到式(1.1)中, 并整理可得

其中

忽略(3.4)式右端的二階小量, 并記ui的近似值為Ui, 可得差分格式HDS1

該格式即為我們熟知的二階中心差分格式(CDS).

Ⅱ 四階橫向差分格式HDS2

由Taylor展式及降階式(2.3)有

將等式(3.7)左右的移項(xiàng)合并, 并整理得

其中

將式(3.9)代入式(3.8), 整理得

其中

將式(3.11)代入到微分方程式(1.1)中, 可得

其中

去掉式(3.13)中小量, 并用近似值Ui代替上式的ui, 可得差分格式HDS2

其截?cái)嗾`差為四階.將(3.15)寫成三點(diǎn)緊湊格式如下

其中

Ⅲ 六階橫向差分格式HDS3

由Taylor展式以及降階式(2.3)有

將式(3.21)中項(xiàng)合并，并整理得

其中

將式(3.23)代入式(3.22), 整理得

其中

將式(3.25)代入到微分方程式(1.1)中, 并整理, 得

其中

忽略式(3.27)右端小量, 可得差分格式HDS3

其截?cái)嗾`差為六階.

4.橫向差分格式的收斂性分析

Ⅰ 中心差分格式(CDS)的收斂性分析

引理4.1(極值原理)設(shè)V={Vi|i ∈是上的網(wǎng)格函數(shù), 差分算子由(4.1)給出,A(x)≥0(x ∈[a,b])成立, 如果

則有

整理得

當(dāng)差分方程組(4.1)為齊次右端與齊次邊界條件時(shí), 由引理4.1知其只有零解, 從而差分格式CDS唯一可解.

引理4.2設(shè)V={Vi|i ∈}是網(wǎng)格上任意的網(wǎng)格函數(shù),差分算子由(4.1)給出,如A(x)≥Amin＞0(x ∈[a,b])成立, 有

證取障礙函數(shù)(barrier function)

簡(jiǎn)單計(jì)算知

令

上式兩端作用差分算子, 有

因而有

證明完畢.

定理4.1Ui是差分格式CDS即(3.6)的解,ui是微分方程(1.1)在結(jié)點(diǎn)處的精確值,當(dāng)A(x)≥Amin＞0(x ∈[a,b])成立, 有

其中M1=‖-Af+A2u+ε(2A′u′+uA′′-f′′)‖∞.

證記誤差為

則由(3.4),(3.6)式有誤差方程

且

故定理4.1成立.

針對(duì)微分方程(1.1), 一般A(x)≥0且函數(shù)A(x),f(x)滿足一定的光滑性, 該方程就是可解的, 中心差分格式可以進(jìn)行數(shù)值計(jì)算.下面我們對(duì)一般情形A(x)≥0(x ∈[a,b])(不一定有A(x)≥Amin＞0)時(shí)中心差分格式的誤差進(jìn)行估計(jì).

引理4.3設(shè)V={Vi|i ∈}是網(wǎng)格上任意的網(wǎng)格函數(shù),差分算子由(4.1)給出,當(dāng)A(x)≥0(x ∈[a,b])時(shí), 有

證取障礙函數(shù)

簡(jiǎn)單計(jì)算知

令

上式兩端作用差分算子, 有

因W(x)＞0,x ∈(a,b), 由引理4.1有

證明完畢.

定理4.2Ui是差分格式CDS即(3.6)的解,ui是微分方程(1.1)在結(jié)點(diǎn)處的精確值,當(dāng)A(x)≥0(x ∈[a,b])時(shí), 有

證在本定理?xiàng)l件下, 首先同定理4.1一樣, 知式(4.16)即仍成立, 再結(jié)合引理4.3得故定理4.2成立.

注4.21)關(guān)于定理4.1的相關(guān)注解對(duì)本定理仍然成立; 2) 當(dāng)A(x)≥Amin＞0 (x ∈(a,b))可使用定理4.1的誤差估計(jì), 當(dāng)僅有A(x)≥0(x ∈[a,b])時(shí)可使用定理4.2的誤差估計(jì).

Ⅱ 橫向差分格式HDS2、HDS3的收斂性分析

利用帶微分余項(xiàng)的泰勒公式計(jì)算, 可得

設(shè)差分方程組為

當(dāng)差分格式HDS2的算子滿足引理4.4的條件時(shí), 顯然齊次差分方程組只有零解, 從而差分格式HDS2唯一可解.

定理成立.

式(3.29)給出了HDS3格式, 泰勒展開(kāi)可知

截?cái)嗾`差為六階, 可類似HDS2一樣進(jìn)行誤差分析, 省略.

5.縱向系列修正差分格式

將降階式(2.3)分別代入式(2.1), (2.2)有

將上兩式(5.1)與(5.2)簡(jiǎn)寫成

其中

聯(lián)立式(5.3),(5.4), 消去得到

上式是無(wú)任何誤差形式精確的差分格式, 下面的關(guān)鍵是如何對(duì)Q1(x),Q2(x),S1(x),S2(x),F1(x),F2(x)保留有限項(xiàng)進(jìn)行截?cái)?

bj(x)、cj(x)、pjk(x)由(2.5)式, 利用Mathematica軟件符號(hào)推導(dǎo)容易寫出其具體表達(dá)式,下面列出其部分式如下.

觀察(5.5)和(5.8)式,Q1(x),Q2(x)中含有級(jí)數(shù)項(xiàng)bj(x)(j=1,2,···), 由遞推公式(2.5), 可算出

觀察(5.6)和(5.9)式, 由遞推公式(2.5), 可算出S1(x),S2(x)中含有的無(wú)窮項(xiàng)cj(x)(j=1,2,···)的表達(dá)式如下

觀察(5.7)和(5.10)式, 對(duì)于F1(x),F2(x)中的pjk(j=k+2,k+3,··· , k= 0,1,2,···), 由遞推公式(2.5), 我們歸納如下.

首先，下面組式是F1(x),F2(x)中關(guān)于f(x)的系數(shù)表達(dá)式pj0(x)(j=2,3,···)

其次,f′(x)的系數(shù)表達(dá)式pj1(x)(j=3,4,5,···)

一般地,f(k)(x)的系數(shù)表達(dá)式pjk(x)(j=k+2,k+3,···)為

Ⅰ 縱向差分格式VDS1

將A(x),f(x)看成常數(shù)A,f,bj(x)、cj(x)、pj0(x)分別保留式(5.12)、(5.13)與(5.14)的第一列作為其近似, 代入(5.5)-(5.10), 整理可得縱向差分格式VDS1

其中i=1,2,··· ,N -1, 且U0=α, UN=β,

上三式代入(5.17), 整理可得VDS1的三點(diǎn)式

Ⅱ 縱向差分格式VDS2

在差分格式VDS1基:上進(jìn)一步修正, 分別保留式(5.12)、(5.13)與(5.14)的前兩列作為bj(x)、cj(x)、pj0(x)的近似, 代入(5.5)-(5.10), 整理可得縱向差分格式VDS2

其中i=1,2,··· ,N -1, 且U0=α, UN=β,

Ⅲ 縱向差分格式VDS3

在差分格式VDS2基:上進(jìn)一步修正, 其中Q1(x),Q2(x),S1(x),S2(x)仍以VDS2中的式(5.22)-(5.25)近似,而(5.7)式F1(x)與(5.10)式F2(x)保留到f′′(x)項(xiàng).對(duì)于pj0(j=2,3,···),pj1(j=3,4,···), 和pj2(j=4,5,···)分別保留第一列, 有

將式(5.22)-(5.25)和式(5.28)-(5.29)代入式(5.11), 可以得到如下差分格式VDS3

6.數(shù)值實(shí)驗(yàn)

例6.1含源反應(yīng)擴(kuò)散方程

精確解為

由圖6.1可見(jiàn), 當(dāng)參數(shù)較小時(shí)該解出現(xiàn)雙邊界層.

圖6.1 參數(shù)不同時(shí)算例6.1的精確解圖

為把舍入誤差影響降低到最小, 以便能考察截?cái)嗾`差實(shí)際精度, 以下各種差分格式結(jié)果均是Mathematica軟件設(shè)定高精度(如ε=0.1精度設(shè)置為50)后編程.

表6.1 算例6.1, 當(dāng)ε=0.1 各種格式的最大誤差

針對(duì)例題1, 因A(x) = 1為常函數(shù), 此時(shí)VDS1與VDS2格式是相同的.由表6.1可看出, 數(shù)值結(jié)果精度由高到低為: HDS3, FD1, HDS2, VDS3, FD2, CDS, VDS1 (VDS2),且HDS3和HDS2數(shù)值結(jié)果明顯優(yōu)于縱向系列差分格式, 說(shuō)明當(dāng)參數(shù)ε較大時(shí), 橫向差分格式HDS3、FD1、HDS2更適用.

減小參數(shù)ε的值, 由表6.2可看出, 數(shù)值結(jié)果精度由高到低約為: VDS3, HDS3, FD1, HDS2,FD2, VDS1 (VDS2), CDS, 縱向格式VDS3的數(shù)值精度優(yōu)于HDS3.

繼續(xù)減小參數(shù)ε的值, 由表6.3可以得到, 數(shù)值方法精度由高到低約為: VDS3, VDS1 (VDS2), HDS3, FD1, HDS2, FD2, CDS, 縱向系列差分格式效果整體優(yōu)于橫向系列差分格式, 說(shuō)明縱向系列差分格式更適用于小參數(shù)反應(yīng)擴(kuò)散問(wèn)題.

表6.2 算例6.1, 當(dāng)ε=0.0001 各種格式的最大誤差

表6.3 算例6.1, 當(dāng)ε=10-7 各種格式的最大誤差

下面給出反應(yīng)項(xiàng)系數(shù)為變系數(shù)的算例實(shí)驗(yàn).

例6.2含源反應(yīng)擴(kuò)散方程

精確解

當(dāng)參數(shù)較小時(shí)該解出現(xiàn)右邊界層.當(dāng)參數(shù)不太小時(shí), 數(shù)值結(jié)果能得到與算例1類似的結(jié)論,篇幅所限, 下面僅給出ε=10-7的實(shí)驗(yàn)結(jié)果.

表6.4 算例6.2, 當(dāng)ε=10-7 各種格式的最大誤差

由表6.4可以得到, 當(dāng)ε=10-7時(shí), 數(shù)值方法精度由高到低為: VDS3, VDS2, VDS1, HDS3,HDS2, CDS, FD1, FD2.因此變系數(shù)情形下仍得到同樣的結(jié)論: 縱向系列差分格式更適宜于求解小參數(shù)反應(yīng)擴(kuò)散問(wèn)題.

由兩個(gè)數(shù)值算例可以看出: 針對(duì)含源反應(yīng)擴(kuò)散問(wèn)題, 中心差分格式雖然不會(huì)出現(xiàn)振蕩現(xiàn)象,但是數(shù)值精度較低.當(dāng)參數(shù)ε不是非常小時(shí), 本文設(shè)計(jì)的橫向系列修正格式HDS3, HDS2適用,但是隨著參數(shù)ε進(jìn)一步減小, 縱向系列差分格式逐漸優(yōu)于橫向系列修正格式, 更適用于小參數(shù)問(wèn)題的求解.

本文設(shè)計(jì)的差分格式可結(jié)合非等距網(wǎng)格剖分, 從而能更細(xì)致地刻畫解在邊界層內(nèi)的變化.另外, 一維降階思想如何應(yīng)用于二維奇異擾動(dòng)問(wèn)題數(shù)值方法設(shè)計(jì), 是我們進(jìn)一步研究的方向.