常系數擴散方程三層十五點差分格式的穩定性

劉瑩 畢卉

摘? 要：基于向后差分格式和Crank-Nicolson格式對二維擴散方程提出一種三層十五點隱式差分格式。采用泰勒展開求出截斷誤差，證明了該格式的相容性，接著用傅里葉變換和Von Neumann條件證明了該格式的無條件穩定性。由于三層差分格式需要兩層啟動條件，在數值實驗中，利用二維Saul′ev差分格式作為三層十五點隱式差分格式的啟動格式。數值試驗表明Saul′ev格式與三層十五點差分格式相結合誤差小，精度高，并且網比的變化對誤差的影響不大。

關鍵詞：三層十五點差分格式；二維擴散方程；穩定性；誤差估計；Saul′ev格式

DOI：10.15938/j.jhust.2023.05.018

中圖分類號： O241.3

文獻標志碼： 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

基金項目： 黑龍江省自然科學基金聯合引導項目（LH2020A015）.

作者簡介：

劉? 瑩（1996—），女，碩士研究生.

通信作者：

畢? 卉（1982—），女，博士，教授，博士研究生導師，E-mail：bihui@hrbust.edu.cn.

0? 引? 言

擴散方程是一類反映液體滲透、氣體擴散、半導體材料雜質擴散等現象的數學模型，在化學、生物、物理等領域都有著非常多的應用。因此，擴散方程數值解法的研究具有重要科學價值。由于有限差分方法（finite difference method，簡稱FDM）簡單靈活且具有很強的通用性，該方法成為一種求解擴散方程數值解的熱門方法［1－9］，越來越多的數值解法［10－11］也被廣泛關注。

基本的差分格式大多是一層和二層的，由于多層差分格式一般具有較高的精確度，本文旨在研究三層差分格式［12－17］。Zhan等［18］用待定系數法構造了求解一維熱傳導方程的三層隱式格式，并研究了該方法的截斷誤差和穩定性條件。Amina等［19］針對擴散方程提出了一種三層九點隱式差分格式，并證明該差分格式與擴散方程相容，二階收斂且無條件穩定。本文基于三層九點差分格式提出一種三層十五點差分格式，并證明該差分格式與二維常系數擴散方程相容且無條件穩定。

論文第1節給出了二維擴散方程的三層十五點差分格式；第2節計算了格式的截斷誤差并判斷了相容性；第3節證明了差分格式的穩定性；第4節通過數值實驗驗證了理論結果；最后在第5節給出了結論。

1? 三層十五點差分格式

二維常系數擴散方程［20］

ut=a2ux2+2uy2，x，y∈R，t>0（1）

其中a是一個正常數。初始條件為

u（x，y，0）=g（x，y），x，y∈R

首先，給出一個近似于微分方程（1）的三層十五點差分格式

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? 截斷誤差

對式（2）中的各項進行泰勒展開，得到

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+…

假設方程（1）的解是光滑的，則有

T（x，y，τ）=13！3ut3njk－a32t22ux2njk+

2uy2njk+h2124ux4njk+h2124uy4njkτ2－

a124ux4njk+a124uy4njkh2+…

因此差分格式（2）的截斷誤差具有二階精度ο（τ2+h2）。由τ→0，h→0得到T（x，y，τ）→0，因此差分格式（2）與微分方程（1）相容。

定理1? 設u（x，t）是定解問題的解，unj是差分格式的解，如果當時間步長τ和空間步長h都趨于零時有

enj=u（xj，tn）－unj→0

那么差分格式是收斂的。

3? 穩定性

把差分格式（2）寫成方便計算的形式

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）

為了證明穩定性，給出了等價于式（10）的兩層方程：

1－23aλ（δ2x+δ2y）un+1jk=

1+23aλ（δ2x+δ2y）vnjk+23aλ（δ2x+δ2y）unjk

vn+1jk=unjk（11）

將δ2xujk，δ2yujk（此處省略上標）代入上述兩層差分方程（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）的特征值函數可以被寫成

μ2－α1－αμ－1+α1－α=0（14）

為了得到格式的穩定性，需要如下引理。

引理1［20］? 實系數二次方程aμ2－bμ－c=0的根按模小于等于1的充要條件是：|b|≤1－c≤2。

已知α≤0，根據式（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）無條件穩定。

定理2? （Von Neumann條件）差分格式穩定的必要條件是當τ≤τ0，nτ≤T，對所有 k∈R有|λj（G（τ，k））|≤1+Mτ，j=1，2，…，p，其中λj（G（τ，k））表示G（τ，k）的特征值，M為常數。

4? 數值算例

已知擴散方程 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，…），其中τ為時間步長，網格比λ=τ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，F1，F2，F3如下：

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作為系數矩陣，于是有

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］計算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）結合來提高精確度。例如，u⌒n+1jk和u⌒n+1jk同時滿足式（15）給定的初邊值條件，則第一層可以由它們的算數平均數U1=12（u⌒1+u⌒1）給出。Saul′ev格式的截斷誤差為Oτh+τ+h2。這里取θ1=θ2=12來計算U1的值，容易證明矩陣A可逆，于是回到式（16）可得到n=2，3，…，N各時間層的數值解。

接下來給出數值解在不同時刻的圖像。圖1和圖2分別給出了當λ=1時在T=1以及T=5時刻對應的數值解。表1給出了在不同網比下不同時刻的數值解與真解。

5? 結? 論

本文討論了常系數擴散方程三層十五點差分格式的穩定性和誤差估計問題。然后用Saul′ev格式求出第一層，再結合三層十五點差分格式求出數值解。

數值結果表明，Saul′ev格式與三層十五點差分格式相結合誤差小，精度高。λ的變化對誤差的影響不大。算法的全部處理表明本文所討論的三層十五點差分方法是無條件穩定、可行的。

參 考 文 獻：

［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］? 孫志忠.偏微分方程數值解法［M］.北京：科學出版社，2005.

［3］? 武莉莉.不可壓磁流體力學方程組的高精度緊致有限差分方法［D］.寧夏：寧夏大學，2021.

［4］? 余德浩， 湯華中.微分方程數值解法［M］.北京：科學出版社，2003.

［5］? 張魯明， 常謙順.復Schrdinger場和實Klein-Gordon場相互作用下一類方程組守恒差分格式的收斂性和穩定性［J］.高等學校計算數學學報，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］? 李小綱.流體力學中雙曲守恒律方程的高精度差分方法研究［D］.西安：西安理工大學，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］? 楊彩杰， 孫同軍.拋物型最優控制問題的Crank-Nicolson差分方法［J］.山東大學學報：理學版，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］.高等應用數學學報，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］.哈爾濱理工大學學報，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］.哈爾濱理工大學學報，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］.山東師范大學學報：自然科學版，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］.四川大學學報：自然科學版，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］.陜西科技大學學報，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方程的三層加權平均線性差分格式［J］.四川大學學報：自然科學版，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］謝建強.一維粘性波動方程的三層緊致差分格式［J］.南昌航空大學學報：自然科學版，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］.四川大學學報：自然科學版，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］詹涌強， 凌婷.求解一維熱傳導方程的一族三層隱格式［J］.西南師范大學學報：自然科學版，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］阿米娜·沙比爾， 楊慶之.常系數擴散方程的三層九點差分格式的穩定性（英文）［J］.高等學校計算數學學報，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］李榮華.偏微分方程數值解法［M］.北京：高等教育出版社，2010.

［21］孫潔.關于向前-向后熱方程的數值方法［D］.杭州：浙江大學，2008.

（編輯：溫澤宇）