劉瑩 畢卉



摘? 要:基于向后差分格式和Crank-Nicolson格式對二維擴(kuò)散方程提出一種三層十五點(diǎn)隱式差分格式。采用泰勒展開求出截?cái)嗾`差,證明了該格式的相容性,接著用傅里葉變換和Von Neumann條件證明了該格式的無條件穩(wěn)定性。由于三層差分格式需要兩層啟動條件,在數(shù)值實(shí)驗(yàn)中,利用二維Saul′ev差分格式作為三層十五點(diǎn)隱式差分格式的啟動格式。數(shù)值試驗(yàn)表明Saul′ev格式與三層十五點(diǎn)差分格式相結(jié)合誤差小,精度高,并且網(wǎng)比的變化對誤差的影響不大。
關(guān)鍵詞:三層十五點(diǎn)差分格式;二維擴(kuò)散方程;穩(wěn)定性;誤差估計(jì);Saul′ev格式
DOI:10.15938/j.jhust.2023.05.018
中圖分類號: O241.3
文獻(xiàn)標(biāo)志碼: A
文章編號: 1007-2683(2023)05-0143-07
Stability of Three-level Fifteen-point Difference Scheme
for Diffusion Equations with Constant Coefficients
LIU Ying,? BI Hui
(School of Sciences, Harbin University of Science and Technology, Harbin 150080, China)
Abstract:Based on backward difference and Crank-Nicolson scheme, a three-level fifteen-point implicit difference scheme for two-dimensional diffusion equation is proposed. The truncation error is obtained by Taylor expansion, and the compatibility of the scheme is proved, then the unconditional stability of the scheme is proved by Fourier transform and Von Neumann condition. Since the three-level difference scheme needs two-level starting conditions, two-dimensional Saul′ev difference scheme is used as the starting scheme of three-level fifteen-point implicit difference scheme in numerical experiments. Numerical experiments show that the combination of Saul′ev scheme and three-level fifteen-point difference scheme has small error and high accuracy, and the change of network ratio has little effect on the error.
Keywords:three-level fifteen-point difference scheme; two-dimensional diffusion equation; stability; error estimation; Saul′ev scheme
收稿日期: 2022-05-22
基金項(xiàng)目: 黑龍江省自然科學(xué)基金聯(lián)合引導(dǎo)項(xiàng)目(LH2020A015).
作者簡介:
劉? 瑩(1996—),女,碩士研究生.
通信作者:
畢? 卉(1982—),女,博士,教授,博士研究生導(dǎo)師,E-mail:bihui@hrbust.edu.cn.
0? 引? 言
擴(kuò)散方程是一類反映液體滲透、氣體擴(kuò)散、半導(dǎo)體材料雜質(zhì)擴(kuò)散等現(xiàn)象的數(shù)學(xué)模型,在化學(xué)、生物、物理等領(lǐng)域都有著非常多的應(yīng)用。因此,擴(kuò)散方程數(shù)值解法的研究具有重要科學(xué)價(jià)值。由于有限差分方法(finite difference method,簡稱FDM)簡單靈活且具有很強(qiáng)的通用性,該方法成為一種求解擴(kuò)散方程數(shù)值解的熱門方法[1-9],越來越多的數(shù)值解法[10-11]也被廣泛關(guān)注。
基本的差分格式大多是一層和二層的,由于多層差分格式一般具有較高的精確度,本文旨在研究三層差分格式[12-17]。Zhan等[18]用待定系數(shù)法構(gòu)造了求解一維熱傳導(dǎo)方程的三層隱式格式,并研究了該方法的截?cái)嗾`差和穩(wěn)定性條件。Amina等[19]針對擴(kuò)散方程提出了一種三層九點(diǎn)隱式差分格式,并證明該差分格式與擴(kuò)散方程相容,二階收斂且無條件穩(wěn)定。本文基于三層九點(diǎn)差分格式提出一種三層十五點(diǎn)差分格式,并證明該差分格式與二維常系數(shù)擴(kuò)散方程相容且無條件穩(wěn)定。
論文第1節(jié)給出了二維擴(kuò)散方程的三層十五點(diǎn)差分格式;第2節(jié)計(jì)算了格式的截?cái)嗾`差并判斷了相容性;第3節(jié)證明了差分格式的穩(wěn)定性;第4節(jié)通過數(shù)值實(shí)驗(yàn)驗(yàn)證了理論結(jié)果;最后在第5節(jié)給出了結(jié)論。
1? 三層十五點(diǎn)差分格式
二維常系數(shù)擴(kuò)散方程[20]
ut=a2ux2+2uy2,x,y∈R,t>0(1)
其中a是一個正常數(shù)。初始條件為
u(x,y,0)=g(x,y),x,y∈R
首先,給出一個近似于微分方程(1)的三層十五點(diǎn)差分格式
un+1jk-un-1jk2τ-a3h2(δ2xun+1jk+δ2xunjk+δ2xun-1jk+δ2yun+1jk+δ2yunjk+δ2yun-1jk)=0(2)
其中:τ為時間步長;h為空間步長;x和y的步長相同,都為h。
δ2xujk=uj+1,k-2ujk+uj-1,k
δ2yujk=uj,k+1-2ujk+uj,k-1
2? 截?cái)嗾`差
對式(2)中的各項(xiàng)進(jìn)行泰勒展開,得到
un+1jk-un-1jk2τ=utnjk+13!3ut3njkτ2+…(3)
δ2xun+1jk=2ux2n+1jkh2+1124ux4n+1jkh4+…=
2ux2njkh2+1124ux4njkh4+t2ux2njk+
h2124ux4njkτh2+12!2t22ux2njk+
h2124ux4njkτ2h2+…(4)
δ2xunjk=2ux2njkh2+1124ux4njkh4+…(5)
δ2xun-1jk=2ux2njkh2+1124ux4njkh4-
t2ux2njk+h2124ux4njkτh2+
12!2t22ux2njk+h2124ux4njkτ2h2+…(6)
δ2yun+1jk=2uy2n+1jkh2+1124uy4n+1jkh4+…=
2uy2njkh2+1124uy4njkh4+t2uy2njk+
h2124uy4njkτh2+12!2t22ux2njk+
h2124ux4njkτ2h2+…(7)
δ2yunjk=2uy2njkh2+1124uy4njkh4+…(8)
δ2yun-1jk=2uy2njkh2+1124uy4njkh4-
t2uy2njk+h2124uy4njkτh2+
12!2t22ux2njk+h2124ux4njkτ2h2+…(9)
把式(3)~式(9)代入差分格式(2),有
un+1jk-un-1jk2τ-a3h2[δ2xun+1jk+δ2xunjk+δ2xun-1jk+
δ2yun+1jk+δ2yunjk+δ2yun-1jk]=utnjk-
a2ux2njk-a2uy2njk+13!3ut3njkτ2-
a124ux4njkh2-a124uy4njkh2-a32t22ux2njk+
2uy2njk+h2124ux4njk+h2124uy4njkτ2+…
假設(shè)方程(1)的解是光滑的,則有
T(x,y,τ)=13!3ut3njk-a32t22ux2njk+
2uy2njk+h2124ux4njk+h2124uy4njkτ2-
a124ux4njk+a124uy4njkh2+…
因此差分格式(2)的截?cái)嗾`差具有二階精度ο(τ2+h2)。由τ→0,h→0得到T(x,y,τ)→0,因此差分格式(2)與微分方程(1)相容。
定理1? 設(shè)u(x,t)是定解問題的解,unj是差分格式的解,如果當(dāng)時間步長τ和空間步長h都趨于零時有
enj=u(xj,tn)-unj→0
那么差分格式是收斂的。
3? 穩(wěn)定性
把差分格式(2)寫成方便計(jì)算的形式
12(un+1jk-un-1jk)=13aλ(δ2xun+1jk+δ2xunjk+δ2xun-1jk+δ2yun+1jk+δ2yunjk+δ2yun-1jk)
其中λ=τh2,整理得到
1-23aλδ2x-23aλδ2yun+1jk=1+23aλδ2x+
23aλδ2yun-1jk+23aλ(δ2x+δ2y)unjk(10)
為了證明穩(wěn)定性,給出了等價(jià)于式(10)的兩層方程:
1-23aλ(δ2x+δ2y)un+1jk=
1+23aλ(δ2x+δ2y)vnjk+23aλ(δ2x+δ2y)unjk
vn+1jk=unjk(11)
將δ2xujk,δ2yujk(此處省略上標(biāo))代入上述兩層差分方程(11),得到
un+1jk-23aλ(δ2xun+1jk+δ2yun+1jk)=vnjk+
23aλ(δ2xvnjk+δ2yvnjk)+23aλ(δ2xunjk+δ2yunjk)
vn+1jk=unjk
如果取Unjk=(unjk,vnjk)T,其中Unjk是一個兩行一列的矩陣,然后把上面的方程寫成向量形式,得到
D1Un+1j+1k+D2Un+1jk+D1Un+1j-1k+D1Un+1jk+1+
D1Un+1jk-1=D3Unj+1k+D4Unjk+D3Unj-1k+
D3Unjk+1+D3Unjk-1(12)
其中
D1=-23aλ0
00? D2=1+83aλ001
D3=23aλ23aλ
00? D4=-83aλ1-83aλ
10
如果Unjk=Vne(i(fjh+gkh)),其中Vn與Un形式相同,都是兩行一列矩陣,那么由式(12)有
D1Vn+1eih(f(j+1)+gk)+D2Vn+1eih(fj+gk)+
D1Vn+1eih(f(j-1)+gk)+D1Vn+1eih(fj+(k+1)g)+
D1Vn+1eih(fj+(k-1)g)=D3Vneih(f(j+1)+gk)+
D4Vneih(fj+gk)+D3Vneih(f(j-1)+gk)+
D3Vneih(fj+(k+1)g)+D3Vneih(fj+(k-1)g)(13)
為了簡化矩陣,令
γf=eihf+e-ihf,γg=eihg+e-ihg,
ω=cos(hf)+cos(hg)
將式(13)消去公因式得到
1-φ0
01Vn+1=φ1+φ10Vn
其中φ=23aλγf+23aλγg-83aλ。
由于
eihf=cos(hf)+isin(hf)
e-ihf=cos(hf)-isin(hf),
于是有
-43aλω+1+83aλ0
01Vn+1=
43aλω-83aλ43aλω+1-83aλ
10Vn
令
α=43aλω-83aλ=
43aλ(cos(hf)+cos(hg))-83aλ
顯然α≤0,得到
1-α001Vn+1=α1+α10Vn
于是得到增長因子
G(τ,k)=1-α001-1α1+α10=
α1-α1+α1-α
10
G(τ,k)的特征值函數(shù)可以被寫成
μ2-α1-αμ-1+α1-α=0(14)
為了得到格式的穩(wěn)定性,需要如下引理。
引理1[20]? 實(shí)系數(shù)二次方程aμ2-bμ-c=0的根按模小于等于1的充要條件是:|b|≤1-c≤2。
已知α≤0,根據(jù)式(14),有b=α1-α,c=1+α1-α
并且|b|≤1-c=1-1+α1-α=-2α1-α且1-c=-2α1-α≤2。由引理1,可以得到|μi|≤1(i=1,2),所以|G|≤1。這樣就滿足了Von Neumann條件,因此差分格式(2)無條件穩(wěn)定。
定理2? (Von Neumann條件)差分格式穩(wěn)定的必要條件是當(dāng)τ≤τ0,nτ≤T,對所有 k∈R有|λj(G(τ,k))|≤1+Mτ,j=1,2,…,p,其中λj(G(τ,k))表示G(τ,k)的特征值,M為常數(shù)。
4? 數(shù)值算例
已知擴(kuò)散方程 ut=4-2(uxx+uyy)
(x,y)∈G=(0,1)×(0,1),t>0
u(0,y,t)=u(1,y,t)=0,0≤y≤1,t>0
u(x,0,t)=u(x,1,t)=0,0≤x≤1,t>0
u(x,y,0)=sinπxsinπy(15)
通過變量分離可以得到方程的解析解
u=sinπxsinπyexp-π28t
(x,y)∈G=(0,1)×(0,1),t>0。
離散方程(15)的定義域:
令xj=jh,yk=kh(j,k=0,1,…,J),tn=nτ(n=1,2,…),其中τ為時間步長,網(wǎng)格比λ=τh2,重新排列(2),得到
un+1jk-23aλun+1j+1k+43aλun+1jk-23aλun+1j-1k-
23aλun+1jk+1+43aλun+1jk-23aλun+1jk-1=un-1jk+
23aλun-1j+1k-43aλun-1jk+23aλun-1j-1k+
23aλun-1jk+1-43aλun-1jk+23aλun-1jk-1+
23aλunj+1k-43aλunjk+23aλunj-1k+
23aλunjk+1-43aλunjk+23aλunjk-1
令j=1:J-1,k=1:J-1
Un=[un11,un21,…,un(J-1)1,un12,un22,…,
un(J-1)2,…,un1(J-1),un2(J-1),…,un(J-1)(J-1)]T
取J-1階方陣Aii,Bii,Cii,F(xiàn)1,F(xiàn)2,F(xiàn)3如下:
Aii=
1+83aλ-23aλ
-23aλ1+83aλ-23aλ
-23aλ1+83aλ-23aλ
-23aλ1+83aλ
Bii=
-83aλ23aλ
23aλ-83aλ23aλ
23aλ-83aλ23aλ
23aλ-83aλ
Cii=
1-83aλ23aλ
23aλ1-83aλ23aλ
23aλ1-83aλ23aλ
23aλ1-83aλ
F1=diag-23aλ, …, -23aλ
F2=diag23aλ, …, 23aλ
F3=diag23aλ, …, 23aλ
記
A=A11F1
F1A22F1
F1AJ-2,J-2F1
F1AJ-1,J-1
B=B11F2
F2B22F2
F2BJ-2,J-2F2
F2BJ-1,J-1
C=C11F3
F3C22F3
F3CJ-2,J-2F3
F3CJ-1,J-1
顯然A,B,C均為(J-1)2階方陣。
用A,B,C作為系數(shù)矩陣,于是有
AUn+1+M=BUn+CUn-1+Q
其中M=m1m2mJ-2mJ-1
Q=q1q2qJ-2qJ-1
m1=-23aλun+10,1-23aλun+11,0
-23aλun+12,0
-23aλun+13,0
-23aλun+1J-2,0
-23aλun+1J,1-23aλun+1J-1,0
q1=23aλun01+un10+23aλun-101+un-110
23aλun20+23aλun-120
23aλun30+23aλun-130
23aλunJ-2,0+23aλun-1J-2,0
23aλunJ,1+unJ-1,0+23aλun-1J,1+un-1J-1,0
m2=-23aλun+10,2000-23aλun+1J,2
q2=23aλun0,2+23aλun-10,200023aλunJ,2+23aλun-1J,2
m3=-23aλun+10,3000-23aλun+1J,3
q3=23aλun0,3+23aλun-10,300023aλunJ,3+23aλun-1J,3
m4=-23aλun+10,4000-23aλun+1J,4
q4=23aλun0,4+23aλun-10,400023aλunJ,4+23aλun-1J,4
一直到mJ-2和qJ-2。
mJ-1=-23aλun+10,J-1+un+11,J-23aλun+12,J-23aλun+13,J-23aλun+1J-2,J-23aλun+1J,J-1+un+1J-1,J
qJ-1=23aλ(un0,J-1+un1,J)+23aλ(un-10,J-1+un-11,J)23aλun2,J+23aλun-12,J23aλunJ-2,J+23aλun-1J-2,J23aλ(unJ,J-1+unJ-1,J)+23aλ(un-1J,J-1+un-1J-1,J)
顯然M=0,Q=0。這樣有
Un+1=A-1BUn+A-1CUn-1(16)
對于三層格式,U0是已知的初始條件,U1未知,這里選擇二維Saul′ev差分格式[21]計(jì)算U1:
u⌒n+1jk=aλ1+(θ1+θ2)aλθ1un+1j+1,k+θ2un+1j,k+1+
(1-θ1)unj+1,k+(1-θ2)unj,k+1+unj-1,k+
unj,k-1-(4-θ1-θ2-1aλ)unj,k(17)
u⌒n+1jk=aλ1+(θ1+θ2)aλθ1un+1j-1,k+θ2un+1j,k-1+
(1-θ1)unj-1,k+(1-θ2)unj,k-1+unj+1,k+
unj,k+1-(4-θ1-θ2-1aλ)unj,k(18)
可以把式(17)、(18)結(jié)合來提高精確度。例如,u⌒n+1jk和u⌒n+1jk同時滿足式(15)給定的初邊值條件,則第一層可以由它們的算數(shù)平均數(shù)U1=12(u⌒1+u⌒1)給出。Saul′ev格式的截?cái)嗾`差為Oτh+τ+h2。這里取θ1=θ2=12來計(jì)算U1的值,容易證明矩陣A可逆,于是回到式(16)可得到n=2,3,…,N各時間層的數(shù)值解。
接下來給出數(shù)值解在不同時刻的圖像。圖1和圖2分別給出了當(dāng)λ=1時在T=1以及T=5時刻對應(yīng)的數(shù)值解。表1給出了在不同網(wǎng)比下不同時刻的數(shù)值解與真解。
5? 結(jié)? 論
本文討論了常系數(shù)擴(kuò)散方程三層十五點(diǎn)差分格式的穩(wěn)定性和誤差估計(jì)問題。然后用Saul′ev格式求出第一層,再結(jié)合三層十五點(diǎn)差分格式求出數(shù)值解。
數(shù)值結(jié)果表明,Saul′ev格式與三層十五點(diǎn)差分格式相結(jié)合誤差小,精度高。λ的變化對誤差的影響不大。算法的全部處理表明本文所討論的三層十五點(diǎn)差分方法是無條件穩(wěn)定、可行的。
參 考 文 獻(xiàn):
[1]? WANG Tingchun, GUO Boling. Analysis of Some Finite Difference Schemes for Two-Dimensional Ginzburg-Landau Equation [J]. Numer Methods Partial Differential Equations, 2011, 27(5): 1340.
[2]? 孫志忠.偏微分方程數(shù)值解法[M].北京:科學(xué)出版社,2005.
[3]? 武莉莉.不可壓磁流體力學(xué)方程組的高精度緊致有限差分方法[D].寧夏:寧夏大學(xué),2021.
[4]? 余德浩, 湯華中.微分方程數(shù)值解法[M].北京:科學(xué)出版社,2003.
[5]? 張魯明, 常謙順.復(fù)Schrdinger場和實(shí)Klein-Gordon場相互作用下一類方程組守恒差分格式的收斂性和穩(wěn)定性[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),2000(4):362.
ZHANG Luming, CHANG Qianshun.Convergence and Stability of a Conservative finite Difference Scheme for a Class of Equation system in Interaction of Complex Schrdinger Field and Real Klein-Gordon Field [J].Journal of Computational Mathematics,2000(4): 362.
[6]? 李小綱.流體力學(xué)中雙曲守恒律方程的高精度差分方法研究[D].西安:西安理工大學(xué),2020.
[7]? ZHANG Luming. Convergence of a Conservative Difference Schemes for a Class of Klein-Gordon-Schrdinger Equations in one Space Dimension [J]. Applied Mathematics and Computation,2000,163(1):343.
[8]? 楊彩杰, 孫同軍.拋物型最優(yōu)控制問題的Crank-Nicolson差分方法[J].山東大學(xué)學(xué)報(bào):理學(xué)版,2020,55(6):115.
YANG Caijie, SUN Tongjun. Crank-Nicolson Finite Difference Method for Parabolic Optimal Control Problem [J]. Journal of Shandong University: Science Edition, 2020, 55 (6):115.
[9]? 王廷春, 張魯明, 陳芳啟, 等.求解Klein-Gordon-Schrdinger方程組的一個新型守恒差分算法的收斂性分析[J].高等應(yīng)用數(shù)學(xué)學(xué)報(bào),2008(1):41.
WANG Tingchun, ZHANG Luming, CHEN Fangqi, et al. Convergence Analysis of a New Conservative Difference Algorithm for Solving Klein-Gordon-Schrdinge Equations[J]. Applied Mathematics A Journal of Chinese, 2008(1):41.
[10]畢卉, 陳莎莎.四階線性方程局部間斷Galerkin方法的誤差估計(jì)[J].哈爾濱理工大學(xué)學(xué)報(bào),2021,26(4):159.
BI Hui, CHEN Shasha. Error Estimates for Local Discontinuous Galerkin Methods for Linear Fourth-order Equations[J]. Journal of Harbin University of Science and Technology,2021,26(4):159.
[11]畢卉, 孟雄, 孫陽.求解雙曲守恒律方程的高階TVD格式[J].哈爾濱理工大學(xué)學(xué)報(bào),2010,15(3):54.
BI Hui, MENG Xiong, SUN Yang. A Higher Order TVD Scheme for Hyperbolic Conservation Laws[J]. Journal of Harbin University of Science and Technology,2010,15(3):54.
[12]蘇保金, 姜子文.二維擬線性粘性波動方程的三層緊致差分格式[J].山東師范大學(xué)學(xué)報(bào):自然科學(xué)版,2019,34(2):171.
SU Baojin, JIANG Ziwen. A Three Level Compact Difference Scheme For Solving A Two-Dumensional Quasilinear Viscous Wave Equation[J]. Journal of Shandong Normal University: Natural Science Edition,2019,34(2):171.
[13]李佳佳, 張虹, 王希,等.Rosenau-KdV-RLW方程的三層線性化差分格式[J].四川大學(xué)學(xué)報(bào):自然科學(xué)版,2018,55(6):1137.
LI Jiajia, ZHANG Hong, WANG Xi, et al. A Three-level Linearized Difference Scheme for Rosenau-KdV-RLW Equation[J].Journal of Sichuan University: Natural Science Edition,2018,55(6):1137.
[14]常紅, 丁丹平.Camassa-Holm方程的一種三層守恒有限差分格式[J].陜西科技大學(xué)學(xué)報(bào),2017,35(3):186.
CHANG Hong, DING Danping.A Three-level Conservative Finite Difference Scheme for Camassa-Holm Equation[J].Journal of Shaanxi University of Science and Technology,2017,35(3):186.
[15]趙紅偉, 胡兵, 鄭茂波.General Improved KdV方程的三層加權(quán)平均線性差分格式[J].四川大學(xué)學(xué)報(bào):自然科學(xué)版,2017,54(1):12.
ZHAO Hongwei, HU Bing, ZHENG Maobo. Three-level Average Linear Difference Scheme for the General Improved KdV Equation[J].Journal of Sichuan University: Natural Science Edition,2017,54(1):12.
[16]謝建強(qiáng).一維粘性波動方程的三層緊致差分格式[J].南昌航空大學(xué)學(xué)報(bào):自然科學(xué)版,2016,30(2):50.
XIE Jianqiang. A Three level Compact Difference Scheme for Solving A One-dimensional Viscous Wave Equation[J].Journal of Nanchang Hangkong University: Natural Science Edition,2016,30(2):50.
[17]杜瑜, 徐友才, 胡兵.Rosenau-Burgers方程的三層差分格式(英文)[J].四川大學(xué)學(xué)報(bào):自然科學(xué)版,2010,47(1):1.
DU Yu,XU Youcai, HU Bing.The Three Level Finite Difference Scheme for Rosenau-Burgers Equation[J].Journal of Sichuan University: Natural Science Edition,2010,47(1):1.
[18]詹涌強(qiáng), 凌婷.求解一維熱傳導(dǎo)方程的一族三層隱格式[J].西南師范大學(xué)學(xué)報(bào):自然科學(xué)版,2020,45(11):1.
ZHAN Yongqiang, LING Ting.A Class of Three-level Implicit Difference Scheme for Solving One-dimensional Heat Conduction Parabolic Equations[J].Journal of Southwest Normal University: Natural Science Edition,2020,45(11):1.
[19]阿米娜·沙比爾, 楊慶之.常系數(shù)擴(kuò)散方程的三層九點(diǎn)差分格式的穩(wěn)定性(英文)[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),2020,42(2):148.
AMINA Shabel, YANG Qingzhi. Stability of Three-level Nine-point Difference Scheme for Constant Coefficient Diffusion Equations[J].Numerical Mathematics A Journal of Chinese Universities, 2020,42(2):148.
[20]李榮華.偏微分方程數(shù)值解法[M].北京:高等教育出版社,2010.
[21]孫潔.關(guān)于向前-向后熱方程的數(shù)值方法[D].杭州:浙江大學(xué),2008.
(編輯:溫澤宇)