求解抛物型方程的一种高精度紧致差分格式

杨晓佳，魏剑英

(宁夏大学数学与计算机学院, 宁夏 银川 750021)



求解抛物型方程的一种高精度紧致差分格式

杨晓佳，魏剑英

(宁夏大学数学与计算机学院, 宁夏 银川 750021)

摘要：利用四阶Padé逼近公式和扩展的1/3-Simpson公式，构造一种求解一维抛物型方程的高精度紧致隐式差分格式，其截断误差为O(τ4+h4). 然后通过理论分析证明此格式是无条件稳定的，并通过数值实验验证本文中格式的精确性和可靠性.

关键词：抛物型方程；高精度紧致格式；无条件稳定；扩展的1/3-Simpson公式；Padé逼近

0引言

抛物型方程的数值解法一直是科研工作者们关注的热点问题之一，该类方程可描述众多物理现象，如大气、海洋、江河及地下水污染物的扩散过程，因此对此类方程模型进行数值研究具有重要意义. 数值求解该类方程的关键问题是如何在计算工作量尽可能小的前提下，使计算精度尽可能得到提高. 目前，已经发展了很多求解该方程的数值方法[1-9]，其中Crank-Nicolson方法是具有二阶精度的绝对稳定的经典方法. 但是Crank-Nicolson方法计算某些问题时会在边界附近出现数值振荡现象[1]；开依沙尔·热合曼[2]采用半离散方法推导了一个三阶精度的半离散隐格式；Sallam等人[3]采用中心差分格式对空间进行离散，得到关于时间t的半离散微分方程后采用了扩展的辛普森(Simpson)公式，得到了一个O(τ4+h2)的隐式紧致差分格式；Chawia等人[4]采用中心差分格式和广义的梯形公式，得到了一个含参数的无条件稳定的隐式紧致差分格式，其截断误差为O(τ2+h2)；马明书[5]基于泰勒级数展开和待定系数的方法构造了一维抛物型方程精度为O(τ3+h4)的3层3点显格式，其稳定性条件为r<1/2(r为网格比)，即格式是条件稳定的；Chawia等人[6]应用含一族参数的L-稳定的Simpson公式得到了一个无条件稳定的含参数的隐格式，其截断误差为O(τ4+h4)，但此格式要依赖参数，并且文献中并未给出选择适当参数的方法；Sun等[7]针对一维传热方程，采用边界值方法和四阶紧致差分格式分别对时间项和空间项进行离散，得到了一个截断误差为O(τ4+h4)的无条件稳定的隐格式；葛永斌等[8]利用一阶微商和二阶微商的四阶紧致差分逼近公式，推导出了求解一维扩散方程的两种精度分别为O(τ2+h4)和O(τ4+h4)的隐式紧致差分格式，其中O(τ2+h4)格式是两层无条件稳定的，而O(τ4+h4)格式是三层无条件不稳定的；杨国锋等[9]利用加权平均的思想和二阶微商的四阶紧致差分逼近公式推导出了一个O(τ2+h4)的半显格式，通过Fourier分析方法证明了该格式是无条件稳定的；詹涌强[10]提出了求解一维抛物型方程的一族两层六点隐式格式，格式的截断误差为O(τ2+h4)，并通过Fourier方法证明了差分格式当1/2≤θ≤1(θ为参数)时，格式是无条件稳定的，当0≤θ<1/2时，只有r≤1/6(1-2θ)时，格式才是稳定的；Hao[11]采用区域分裂的方法，构造了一种求解一维和二维抛物型方程的高精度显格式，其截断误差为O(τ3+h3)，并且通过理论分析的方法证明其格式是无条件稳定的；文献[12]中对空间变量应用中心差分格式和紧致差分格式离散，时间变量采用二级四阶Runge-Kutta方法，构造了求解一维扩散方程的精度为O(τ4+h2)和O(τ4+h4)的两种无条件稳定的隐式差分格式，但这两种格式都只能计算边界为零的问题；詹涌强等人[13]利用抛物型方程解的一阶偏导数在特殊节点处的一个差分近似式和二阶中心差分近似式，并结合待定系数法构造了一族隐式差分格式，格式的截断误差为O(τ3+h4)，稳定性条件为r>1/6.

本文中利用四阶Padé逼近公式和扩展的1/3-Simpson公式，构造出一种求解一维抛物型方程的高精度隐式紧致差分格式，然后对格式的稳定性进行理论分析. 最后，通过数值算例对本文中方法的精确性和可靠性进行验证.

1差分格式的建立

考虑一维抛物型方程的初边值问题

ut=vuxx,00

(1)

给定初始条件为

u(x,0)=φ(x),0

(2)

给定边界条件为

u(0,t)=a(t),u(l,t)=b(t),t>0

(3)

其中u(x,t)是待求未知函数，v为扩散项系数，a(t)、b(t)、φ(x)均为已知函数.

首先对计算区域进行离散，设时间步长为τ，空间方向N等分网格，空间步长为h=l/N，网格节点为(xj,tn)，其中xj=jh,tn=nτ,j=0,1,2,…,N,n>0，对于任意取定的t，设uj(t)是u(xj,t)的近似值，(uxx)j(t)是uxx(xj,t)的近似值，(ut)j(t)是ut(xj,t)的近似值，依此类推.

对空间内部节点采用四阶精度的Padé公式[14]逼近，则有

(4)

对于空间边界节点的处理，由式(1)与式(3)可得

(5)

(6)

令

(7)

由(7)式可得

(8)

其中

E=J-1S,C(t)=J-1[C2(t)-h2C1(t)].

由(1)式和(8)式可得

(9)

(10)

其中α=n-1,n-1/2,n,则有

(11)

其中

(12)

(13)

将(10)式代入(12)式和(13)式可得

(14)

(15)

整理(14)式和(15)式可得

(16)

(17)

其中I为单位矩阵.

(18)

(19)

将(18)式与(19)式相加可得

(20)

(20)式即为所构造的差分格式. 通过对空间方向的离散过程不难发现，其空间具有四阶精度. 另外，由文献[3]中的研究结果可知，时间方向也具有四阶精度. 因此，本文中所构造的差分格(20)式的截断误差为O(τ4+h4).

2稳定性分析

下面分析格式(20)的稳定性，假定边界条件精确满足.

引理 2.1假设λ和x∈RN-1(x≠0)分别是矩阵J-1S的特征值和特征向量(J和S如上文定义)，则λ是实数且有λ>0.

引理 2.1的证明由于λ和x分别是矩阵J-1S的特征值和特征向量，那么

J-1Sx=λx⟹λxTJx=xTSx.

令

x=(x1,x2,…,xN-1)T,

那么

所以

xTSx>0.

而

则

xTJx>0.

由

λxTJx=xTSx⟹λ>0,

且λ为实数.

引理 2.2的证明当z>0时，有

288z+24z3>0⟹

144z+12z3>-144z-12z3⟹

z4+60z2+144z+12z3+144>z4+60z2-144z-12z3+14⟹4

(12-6z+z2)2<(12+6z+z2)2⟹

(1-z/2+z2/12)2<(1+z/2+z2/12)2,

由此可得

定理2.1格式(20)式是无条件稳定的.

定理2.1的证明假设λi(i=1,2,…)是矩阵E的任意特征值，则

是

的特征值，根据引理2.1，我们可以得到矩阵E的所有特征值都为正的实数，再根据引理2.2可得

因此，格式(20)是无条件稳定的.

3数值算例

为了验证所提格式的精确性和可靠性，现考虑如下3个有精确解的初边值问题.

收敛阶

其中Error(N1)和Error(N2)分别为网格数为N1和N2时对应的最大绝对误差，Uj为数值解，uj为精确解.

算例1[7]

u(x,0)=sin(πx),0

u(0,t)=u(1,t)=0,t>0.

该问题的精确解为

u(x,t)=e-π2tsinπx.

算例2

u(x,0)=sin(x)+x,0

u(0,t)=0,u(π,t)=π,t>0,

该问题的精确解为

u(x,t)=x+e-tsinx.

算例3[3]

u(0,t)=0,u(1,t)=1,t>0,

该问题的精确解为

算例4[4]

ut=uxx,00,

u(x,0)=1,0

u(0,t)=u(2,t)=0,t>0,

该问题的精确解为

表1给出了问题1当τ=h时，在t=1时刻对不同N的最大绝对误差及收敛阶，从表中的计算结果可以看出，文献[3]和文献[13]中空间只有二阶精度，文献[7]和本文中格式具有四阶精度，而本文中格式的最大绝对误差比文献[3,7,13]中的要小几个数量级. 表2给出了问题1当τ=h2时，在t=1时刻对不同N的最大绝对误差及收敛阶，文献[4]中是含参数的格式，在计算所有数值算例时均取参数α=0.2，从表中的计算结果可以看出，文献[3]和文献[4]中空间只有二阶精度，文献[13]中空间具有四阶精度，对比表1和表2可知，文献[13]中的方法当τ=h时，空间只有二阶精度，当τ=h2时，空间才能达到四阶精度，而本文中格式无论τ=h，还是τ=h2，空间都具有四阶精度，且其最大绝对误差比文献[3-4,13]中的小几个数量级. 表3给出了问题1当N=32时，在t=1，网格比r=0.8,3.2,12.8,51.2时的最大绝对误差，从表中计算结果可以看出计算均是稳定的，这表明格式是无条件稳定的，这与本文中的理论分析结果是一致的. 表4给出了问题1当h=0.1,τ=0.1时，在t=1时刻不同节点的绝对误差和最大绝对误差，从表中计算结果可以看出，文献[12]中的格式(8)比本文格式的绝对误差大一个数量级，文献[12]中的格式(9)与本文中格式的绝对误差都在同一个数量级上，但本文格式的绝对误差更小一点，而且文献[12]中的两种格式都只能计算边界为零的问题，本文中格式采用了一种较巧妙的边界处理方法，克服了文献[12]中的缺陷；文献[13] 中的绝对误差比本文格式的绝对误差大两个数量级，这表明文中格式的计算结果更为精确.

表1　问题1当τ=h，在t=1时刻对不同N的最大绝对误差及收敛阶

表2　问题1当τ=h2,在t=1时刻对不同N的最大绝对误差及收敛

表3　问题1当N=32时，在t=1时刻对不同r的最大绝对误差

表4　问题1当h=0.1,τ=0.1时，在t=1时刻的绝对误差

表7给出了问题3在t=1时，取不同的h,τr时的最大绝对误差，对比文献[3-4,13]中的计算结果可以发现，本文中格式的计算结果更加精确. 表8给出了问题4当N=32,τ=h2时，在t=5时刻不同节点处的绝对误差和精确解，从表中的数值结果可以看出，本文中格式比文献[3-4,13]中的计算结果更为精确.

表5　问题2当τ=h，在t=1时刻对不同N的最大绝对误差及收敛阶

表6　问题2当τ=h2，在t=1时刻对不同N的最大绝对误差及收敛阶

表7　问题3在t=1时刻的最大绝对误差

表8　问题4当N=32,τ=h2时，在t=5时刻的绝对误差和精确解

4结论

针对一维抛物型方程提出了一种新的隐式高精度紧致差分格式，格式的截断误差为O(τ4+h4)，并通过理论分析证明了格式是无条件稳定的. 通过数值实验将本文中格式与其他文献中格式的计算结果进行了比较，研究表明本文中格式的计算结果更加精确，精度和稳定性与理论分析相吻合，从而验证了本文中格式的精确性与可靠性. 笔者给出了一维抛物型方程的高精度紧致隐格式，该方法可以直接推广到二维和三维问题，我们将另 撰文报道.

5参考文献

[1] Richtmyer R D, Morton K W. Difference methods for initial-value problems[M]. New York：Interscience. 1967.

[2] 开依沙尔·热合曼. 一种一维扩散方程三阶精度的半离散隐式差分格式[J].甘肃联合大学学报：自然科学版, 2009, 23(1): 33-36.

[3] Sallam S, Naim Anwar M, Abdel-Aziz M R . Unconditionally stableC1- cubic spline collocation method for solving parabolic equations[J]. Intern J Computer Math, 2004, 81(7): 813-821.

[4] Chawla M M, Al-Zanaidi M A, Evans D J. Generalized trapezoidal formulas for parabolic equations[J]. Intern J Computer Math, 1998, 70:429-443.

[5] 马明书. 一维抛物型方程的一个新的高精度显式差分格式[J].数值计算与计算机应用, 2001, 22(2): 156-160.

[6] Chawla M M, Al-Zanaidi M A, Al-Shammeri A Z. High-accuracy finite-difference schemes for the diffusion equation[J]. Neural Parallel Sci Comput, 1998( 6): 523-535.

[7] Sun H W, Zhang J. A high compact boundary value method for solving one dimensional heat equations [J]. Numer Method Partial Differen Eq, 2003, 19: 846-857.

[8] 葛永斌, 田振夫, 詹咏. 求解扩散方程的一种高精度隐式差分方法[J].上海理工大学学报, 2005, 27(2): 107-112.

[9] 杨国锋, 李波, 田巧娴, 等. 求解一维抛物型方程的一种高精度半显式差分方法[J].西北师范大学学报：自然科学版, 2008,44(3): 8-11.

[10] 詹涌强. 求解抛物型方程的一族六点隐式差分格式[J].安徽大学学报：自然科学版, 2012, 36 (4): 26-29.

[11] Hao W R , Zhu S H . Domain decomposition schemes with high-order accuracy and unconditional stability[J]. Appl Math Comput, 2013, 219: 6170-6181.

[12] 开依沙尔·热合曼, 努尔买买提·黑力力. 求解扩散方程的二级四阶隐式Runge-Kutta方法[J].湖北大学学报：自然科学版, 2014, 36(5): 476-480.

[13] 詹涌强, 张传林. 解抛物型方程的一族高精度隐式差分格式[J]. 应用数学和力学, 2014, 35 (7): 790-797.

[14] Lele S K. Compact finite difference schemes with spectral-like resolution[J]. J Comput Phys, 1992, 103: 16-42.

[15] Sallam S, Naim Anwar M. Stabilized cubicC1-spline collocation method for solving first-order ordinary initial value problems[J]. Intern J Computer Math, 2000, 74: 87-96.

(责任编辑赵燕)

A high-order compact difference scheme for the parabolic equation

YANG Xiaojia, WEI Jianying

(College of Mathematics and Computer Science,Ningxia University,Yinchuan 750021,China)

Abstract:A high accurate implicit compact difference scheme is constructed for solving the one dimentional parabolic equation by the fourth-order Padé formula combined with time extension of the 1/3-Simpson formulas. The truncation error of the scheme is O(τ4+h4).The stability of the scheme is analyzed by theoritical method. It is unconditionally stable. Numerical experiments are conducted to verify the accuracy and reliability of the present scheme.

Key words:parabolic equation; high-order compact difference scheme; unconditionally stable; extension of the 1/3-Simpson formulas; Padé approximation

中图分类号:O241.82

文献标志码:A

DOI:10.3969/j.issn.1000-2375.2016.02.013

文章编号:1000-2375(2016)01-0160-08

作者简介:杨晓佳(1988-)，女，硕士生，E-mail:yang_xiaoj@sina.com;魏剑英，通信作者，教授

基金项目:宁夏大学自然科学基金 (ZR1407)资助

收稿日期:2015-06-23