祁应楠

(宁夏师范学院数学与计算机科学学院，宁夏固原 756000)



求解一维抛物型方程的高精度有限差分方法

祁应楠

(宁夏师范学院数学与计算机科学学院，宁夏固原 756000)

针对一维抛物型方程，采用样条函数近似和Padé公式，构造了一种高精度有限差分格式.该格式关于时间和空间均具有六阶精度，并且从理论上被证明是无条件稳定的.通过数值实验验证了本文方法的精确性和稳定性，与文献计算结果比较显示,本文格式的计算结果更加精确.

抛物型方程；样条函数；Padé逼近；高精度；有限差分法

抛物型方程是一类非常重要的偏微分方程，在扩散、热传导、渗流等领域中有着广泛的应用，如绝热过程中气体通过多孔介质的流动问题、一个内部有热源的热传导过程、污染物浓度的扩散问题等都可以用抛物型方程来描述.由于实际问题往往很难得到精确解，因此研究其数值解法具有重要的理论意义和实际应用价值.

对抛物型方程很多情况下都采用有限差分法进行求解.由于传统的差分格式往往精度较低，一般仅具有一阶或二阶精度，因此不少研究者致力于发展高精度差分格式.马明书等[1]针对一维抛物型方程建立了一族隐式差分格式，其时间为三阶精度，空间为六阶精度，格式是无条件稳定的；Sun 和Zhang[2]针对一维抛物型方程，对空间导数项利用四阶紧致差分，对时间项利用边界值方法，构造了无条件稳定的四阶精度差分方法.Sallam等[3]对空间变量利用中心差分公式进行离散，得到关于时间变量的常微分方程组，然后基于C1三次样条配置法得到了时间为四阶精度而空间为二阶精度的无条件稳定的差分格式；葛永斌等[4]针对一维抛物型方程，利用Crank-Nicolson格式的思想和二阶导数的四阶紧致差分公式，分别建立了时间为二阶、空间为四阶和时间空间均为四阶精度的紧致隐式差分格式，并利用离散Fourier分析法证明了前一种差分格式是无条件稳定的，但后一种差分格式是无条件不稳定的.Wen和Shao[5]采用区域分解法，构造了一种求解一维和二维抛物型方程的高精度显格式，其截断误差为O(τ3+h3)，并且是无条件稳定的.Gao和Sun[6]建立了Neumann边界条件下一维抛物型方程的一种紧致差分格式，此格式在内网格点上的计算精度为时间二阶和空间四阶，而在边界点上的计算精度为时间二阶和空间三阶.刘蕤和高锐敏[7]基于广义梯形公式和四次样条函数，建立了Neumann边界条件下的一维抛物型方程的一族含参数的隐式差分格式，该格式在空间方向上达到了四阶精度，但时间方向上仅有二阶精度.当取特殊参数时，该格式的时间精度可提高到三阶.陈国玉等[8]建立了一维抛物型方程的一种条件稳定的三层高精度差分格式，在时间方向上具有二阶精度，而空间方向上具有四阶精度.

本文针对一维抛物型方程，采用样条函数近似和Padé公式，构造一种时间和空间均具有六阶精度的有限差分格式，并且从理论和数值实验两方面分析和验证本文所提格式的精确性和稳定性.考虑一维抛物型方程

(1)

初始条件

(2)

边界条件

(3)

其中φ(x)是已知函数，a>0为扩散项系数，α,β是常数.

1 高阶有限差分格式

(4)

(5)

(6)

(7)

根据文献[9]的结论，有

(8)

(9)

组合(1),(8)和(9)式，即可得到

(10)

即

(11)

又因为

(12)

于是有

(13)

即

(14)

舍掉高阶截断误差项，则方程(1)的半离散格式可写为

(15)

(15)式包含N-3个方程，但根据计算域离散有N-1个未知节点，因此还需补充两个边界点方程，方程组才有唯一解.为了得到与(15)式相匹配的六阶精度格式，定义以下恒等式[10]：

(16)

(17)

利用泰勒级数展开，可以得到以下边界方程：

4u0-7u1+2u2+u3=

(18)

(19)

利用(18)和(19)式，可以得到问题(1)～(3)的两个边界方程分别为：

(20)

(21)

以上两种格式的截断误差均为O(h6).

由(3)式得：u0(t)=α,uN(t)=β,(ut)0(t)=0,(ut)N(t)=0.令

结合(15),(20)和(21)式，有

(22)

其中

由于矩阵A是严格对角占优的，B是不可约对角占优矩阵，那么矩阵A和B都是非奇异的[11]，结合初始条件(2)，有

(23)

其中

则(23)式的解为

(24)

其中

I是(N-1)×(N-1)阶单位矩阵,则(24)式可以变形为

(25)

令X(t)=U(t)-B-1b,X0=U0-B-1b,C=A-1B,则(25)式可重新写为

(26)

令τ为时间步长，则(26)式在点(xi,tn)处的离散结果为

(27)

其中tn=nτ,n=0,1,2,….令X(tn+1)=Xn+1,X(tn)=Xn，则

(28)

其中n=0,1,2,….由于e-τ C是关于矩阵τC的一个无穷级数，为了逼近e-τ C，我们采用如下(3,3)Padé逼近[12]：

(29)

即采用

(120I-60τC+12τ2C2-τ3C3)×

(30)

来逼近e-τ C，那么(28)式的离散格式为

(31)

其中n=0,1,2,….从格式构造过程易知格式(31)的截断误差为O(τ6+h6).

2 格式的稳定性分析

定理1 格式(31)是无条件稳定的.

D(120I-60τC+12τ2C2-τ3C3)×

那么矩阵D的任意特征值为

记

Ej=120-60τ(aj+bji)+12τ2(aj+bji)2-

Fj=120+60τ(aj+bji)+12τ2(aj+bji)2+

则有

于是有

由于aj>0，因此有

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

3 数值算例

为了验证本文格式(31)的性能，我们选用如下几个问题进行数值实验.收敛阶、最大和平均绝对误差的定义如下：

其中，Max_erri(i=1,2)是空间网格步长为hi(i=1,2)时的最大绝对误差，Ave_erri(i=1,2)是空间网格步长是hi(i=1,2)时的平均绝对误差，u(xj)是在xj点处的精确解，Uj是在xj点处的数值解.

问题1[2]ut=(1/π2)uxx,00,u(x,0)=sin(πx),u(0,t)=u(1,t)=0.该问题的精确解为u(x,t)=e-tsin(πx).

表1 问题1在τ=h,t=1.0时的最大绝对误差及收敛阶

表2 问题1在h=1/40,t=1.0时对不同网格比λ的最大绝对误差

表1给出了空间步长h分别取1/8,1/16,1/32,1/64时，问题1在时间t=1.0时刻C-N格式、文献[2]格式和本文格式的最大绝对误差及收敛阶.可以看出，在相同网格等分数下，由本文格式得到的最大绝对误差比C-N格式和文献[2]格式小好几个数量级.当τ=O(h)时，本文格式达到了六阶精度，而C-N格式为二阶精度，文献[2]格式为四阶精度,且只有在τ=O(h2)时精度才能达到四阶.

表2给出了问题1在不同网格比λ下的最大绝对误差，以此来验证三种格式的稳定性.可以看出，计算均是稳定的,并且本文格式的计算结果比C-N格式和文献[2]格式精确得多.另外，从理论上讲本文方法时间和空间均具有六阶精度，由于本文空间导数采用了样条方法(本质上是多项式逼近)，而时间导数采用了(3,3)Padé近似(本质上是指数型逼近)，因此关于时间和空间的截断误差主项系数具有不同的函数特征，从而不能保证计算精度正好是六阶，所以会出现一定的波动现象.

问题2[3]ut=uxx,0

表3 问题2在λ=1/2,t=2.0时的平均绝对误差及收敛阶

表4 问题2在λ=1,t=2.0时的平均绝对误差及收敛阶

表3和表4分别给出了当λ=τ/h2=1/2和1，空间步长h分别取π/8,π/16,π/32,π/64时，问题2在时间t=2.0时刻C-N格式、文献[3]格式和本文格式的平均绝对误差及收敛阶.从表中可以看出，本文格式的结果比C-N格式和文献[3]格式精确得多.本文格式达到了六阶精度，而C-N格式和文献[3]格式仅为二阶精度.

问题3[4]ut=auxx,00,u(x,0)=sin(πx),u(0,t)=u(1,t)=0.该问题的精确解为u(x,t)=e-aπ2tsin(πx).

表5给出了问题3当a=1,τ=h2,h分别取1/8,1/16,1/32,1/64时，在时间t=0.25时刻C-N格式、文献[4]格式和本文格式的最大绝对误差及收敛阶.可以看出，C-N格式仅有二阶精度，文献[4]格式具有四阶精度，而本文格式达到了六阶精度.表6给出了问题3当a=1,h=1/32时不同时间步长τ在t=0.5时刻的最大绝对误差.可以看出，对该问题，网格比λ从0.8变化到51.2，三种格式的计算均是稳定的,并且本文格式的计算结果比C-N格式和文献[4]格式精确得多.

表5 问题3在a=1,τ=h2,t=0.25时的最大绝对误差及收敛阶

表6 问题3在a=1,h=1/32,t=0.5对不同网格比λ的最大绝对误差

4 结论

本文针对一维抛物型方程，首先对空间二阶导数采用六次样条函数进行逼近，从而得到了半离散格式；接下来利用常微分方程精确解公式，得到关于时间t的指数矩阵表达式；然后采用(3,3)Padé逼近得到了一种求解一维抛物型方程的六阶精度有限差分格式，并从理论上证明了该格式是绝对稳定的.最后，数值实验验证了本文格式在时间和空间上的精度都达到了六阶，并且比传统的C-N格式和文献[3-5]格式的计算结果更加精确.

[1] 马明书,李改弟.抛物型方程的一个新的高精度恒稳定的隐式差分格式[J].数学的实践与认识,2001,31(3):355.

[2] UN H W,ZHANG J.A high-order compact boundary value method for solving one-dimensional heat equations[J].NumericalMethodsforPartialDifferentialEquations,2003,19:846.

[3] SALLAM S,NAIM A M,ABDEL-AZIZ M R.Unconditionally stableC1-cubic spline collocation method for solving parabolic equations[J].InternationalJournalofComputerMathematics,2004,81(7):813.

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

[5] WEN R H,SHAO H Z.Domain decomposition schemes with high-order accuracy and unconditional stability[J].AppliedMathematicsandComputation,2013,219:6170.

[6] GAO G H,SUN Z Z.Compact difference schemes for heat equation with Neumann boundary conditions(Ⅱ)[J].NumericalMethodsforPartialDifferentialEquations,2013,29:1459.

[7] 刘蕤,高锐敏.带Neumann边界条件的抛物型方程的样条差分方法[J].郑州大学学报(理学版),2013,45(3):37.

[8] 陈国玉,谢英超,程燕.热传导方程的一个三层差分格式[J].四川兵工学报,2014,35(7):143.

[9] KHAN A,KHAN I,AZIZ T.Sextic spline solution of a singularly perturbed boundary-value problems[J].AppliedMathematicsandComputation,2006,181:432.

[10] RASHIDINIA J,JALILIAN R,MOHAMMADI R,et al.Sextic spline method for the solution of a system of obstacle problems[J].AppliedMathematicsandComputation,2007,190:1669.

[11] 孙文瑜,杜其奎,陈金如.计算方法[M].北京:科学出版社,2007.

[13] 金升平,熊方方,李琼.矩阵的实特征值为正的条件与判断[J].重庆理工大学学报(自然科学版),2015,25(1):117.

(责任编辑 马宇鸿)

A high order finite difference method for solving the one-dimensional parabolic equation

QI Ying-nan

(College of Mathematics and Computer Science,Ningxia Normal University,Guyuan 756000,Ningxia,China)

A high order finite difference scheme for solving the one-dimensional parabolic equation is developed based on the Sextic spline functions for the spatial derivative and (3,3) Padé approximation for the temporal derivative.The scheme is sixth order accuracy in both time and space,and the performance of unconditional stability is proved in theory.The accuracy and effectiveness of the present method are validated by three numerical experiments,the results show that the present method is more accurate than the methods in the literature.

parabolic equation;spline function;Padé approximation;high order accuracy;finite difference method

10.16783/j.cnki.nwnuz.2016.06.006

2015-12-03;修改稿收到日期:2016-06-03

宁夏高等学校科学研究资助项目(NGY2015115);宁夏自然科学基金资助项目(NZ15259，NZ16251);宁夏师范学院科研项目(NXSFZD1603)

祁应楠(1978—)，男，宁夏固原人，副教授，硕士.主要研究方向为偏微分方程数值解法、计算流体力学.

E-mail:421270380@qq.com

O 241.82

A

1001-988Ⅹ(2016)06-0029-06