常系数扩散方程三层十五点差分格式的稳定性

刘莹 毕卉

摘  要：基于向后差分格式和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.

（编辑：温泽宇）