楊曉佳,王燕

(寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院,寧夏 銀川　750021)

?

求解擴(kuò)散方程的高精度顯式緊致差分格式

楊曉佳,王燕

(寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院,寧夏 銀川750021)

首先針對一維擴(kuò)散方程，空間方向采用二階導(dǎo)數(shù)的四階緊致差分公式進(jìn)行離散，時間方向采用泰勒級數(shù)展開的方法進(jìn)行離散，推導(dǎo)出了一種高精度顯式緊致差分格式；然后通過Fourier分析方法給出了格式的穩(wěn)定性條件為λ≤1/2(λ為網(wǎng)格比)；最后通過數(shù)值實(shí)驗(yàn)驗(yàn)證了格式的精確性和可靠性.

擴(kuò)散方程;高精度;緊致格式;顯格式;有限差分法

MSC 2010:35K20

擴(kuò)散方程是一類描述物理量隨時間的擴(kuò)散和衰減規(guī)律的拋物型偏微分方程.自然環(huán)境、工程設(shè)備及生物機(jī)體中的許多物理現(xiàn)象，諸如氣體的擴(kuò)散、液體的滲透等都可用擴(kuò)散方程來描述.由于物理問題本身的復(fù)雜性，其精確解往往不容易求得，因此研究其數(shù)值求解方法無疑具有非常重要的理論意義和實(shí)際應(yīng)用價值.求解擴(kuò)散方程的數(shù)值方法一般可分為顯式方法和隱式方法2種.

近幾年來，研究者們構(gòu)造了許多高精度隱式方法.文獻(xiàn)[1]基于Hermite插值多項(xiàng)式的構(gòu)造思路，推導(dǎo)出了一維含源項(xiàng)擴(kuò)散方程的高精度隱式緊致差分格式；文獻(xiàn)[2]利用一階微商和二階微商的四階緊致差分逼近公式，推導(dǎo)出了求解一維擴(kuò)散方程的2種精度分別為O(τ2+h4)和O(τ4+h4)隱式差分格式，其中O(τ2+h4)格式是2層無條件穩(wěn)定的，而O(τ4+h4)格式是3層無條件不穩(wěn)定的格式.文獻(xiàn)[3]利用拋物型方程解的一階偏導(dǎo)數(shù)在特殊節(jié)點(diǎn)處的一個差分近似式和二階中心差商近似式用待定系數(shù)法構(gòu)造了一族隱式差分格式，格式的截?cái)嗾`差為O(τ3+h4)，穩(wěn)定性條件為r>1/6.文獻(xiàn)[4]用待定系數(shù)法構(gòu)造了求解拋物型方程的一個高精度隱式格式，格式的截?cái)嗾`差達(dá)到O(τ4+h4)，證明了當(dāng)r>1/12時，差分格式是穩(wěn)定的.文獻(xiàn)[5]針對一維拋物型方程的初邊值問題，采用待定系數(shù)法和泰勒級數(shù)展開的方法構(gòu)造了一個2層8點(diǎn)隱式差分格式，其格式的截?cái)嗾`差為O(τ3+h5)，穩(wěn)定性條件是0.001

(1)

初始條件為

(2)

邊界條件為

(3)

其中u(x,t)是待求未知函數(shù)，f(x,t)為非齊次項(xiàng)，a(t)、b(t)、φ(x)均為已知函數(shù).

1　差分格式的建立

考慮方程(1)在第n時刻的值，對空間內(nèi)部節(jié)點(diǎn)采用四階精度的導(dǎo)數(shù)型緊致差分公式[15]，則有

(4)

對于空間邊界節(jié)點(diǎn)的處理，由式(1)與式(3)可得

(5)

(6)

(7)

(8)

由于

(9)

因此

(10)

(11)

(12)

2　穩(wěn)定性分析

下面采用Fourier分析方法對上述建立的差分格式的穩(wěn)定性進(jìn)行分析.假設(shè)邊界條件精確成立，且格式的源項(xiàng)f無誤差存在.

(13)

(14)

對式(13)和式(14)化簡整理得

(15)

(16)

將式(15)代入式(16)得

(17)

從而可得式(12)的增長因子為

(18)

要使式(12)穩(wěn)定，則需|G(τ,ω)|≤1，即

(19)

(20)

式(20)的解為λ∈R.

3　數(shù)值實(shí)驗(yàn)

為了驗(yàn)證本文格式的有效性和可靠性，現(xiàn)考慮如下4個具有精確解的初邊值問題，其中邊界條件給定為Dirichlet邊界條件.采用Fortran 77語言進(jìn)行編程，在Pentium 4個人計(jì)算機(jī)上采用雙精度制進(jìn)行計(jì)算.分別采用文獻(xiàn)[2]、文獻(xiàn)[9]、文獻(xiàn)[10]和本文格式算法進(jìn)行了計(jì)算，給出了不同時刻、不同時間步長、不同空間步長以及不同網(wǎng)格比λ下的最大絕對誤差、CPU時間和收斂階，其中最大絕對誤差及收斂階定義如下：

其中，Error1和Error2分別為空間網(wǎng)格步長為h1和h2時的最大絕對誤差.

算例1[9-10]

問題1是齊次的擴(kuò)散問題，其中右邊界條件是關(guān)于時間t的函數(shù)，從表1計(jì)算結(jié)果可以看出，本文格式對邊界的處理方法是有效的.表1給出了當(dāng)t=1時，分別取網(wǎng)格比λ=1/4和λ=1/2，采用文獻(xiàn)[9]、文獻(xiàn)[10]的方法和本文格式的計(jì)算結(jié)果.文獻(xiàn)[9]中的格式含參數(shù)η，從文中的穩(wěn)定性理論分析可以看出，其格式的穩(wěn)定性條件與參數(shù)有關(guān)，在接下來的計(jì)算中，均取η=-4.當(dāng)網(wǎng)格比λ=1/4時，文獻(xiàn)[9]、文獻(xiàn)[10]的方法和本文格式都具有四階精度，但本文格式的最大絕對誤差小于文獻(xiàn)[9]和文獻(xiàn)[10]中的最大絕對誤差；當(dāng)網(wǎng)格比λ=1/2時，從表中可以看出，當(dāng)網(wǎng)格數(shù)N=16時，文獻(xiàn)[9]的計(jì)算結(jié)果就發(fā)散了；采用文獻(xiàn)[10]所得的最大絕對誤差要比本文格式大3～6個數(shù)量級，而且文獻(xiàn)[10]的方法只具有二階精度，但本文格式仍具有四階精度.表2給出了網(wǎng)格比λ=3/8和λ=1/2時，N=10,n=200(n=t/τ為總的時間步數(shù))時問題1的數(shù)值解和精確解，從表中可以看出本文格式的計(jì)算結(jié)果較文獻(xiàn)[9]和文獻(xiàn)[10]的計(jì)算結(jié)果更精確.同時可以看出，文獻(xiàn)[9]中的結(jié)果發(fā)散.‘div’表示結(jié)果發(fā)散，下同.

表1　問題1在t=1時刻的最大絕對誤差及收斂階

表2　問題1在N=10的數(shù)值解和精確解

算例2

表3給出了問題2在t=1時，分別取網(wǎng)格比λ=1/4和λ=1/2，采用文獻(xiàn)[9]、文獻(xiàn)[10]的方法和本文格式的最大絕對誤差及收斂階，從表中的計(jì)算結(jié)果可以看出，當(dāng)網(wǎng)格比λ=1/4，文獻(xiàn)[9]、文獻(xiàn)[10]的方法和本文格式都具有四階精度，但本文格式的最大絕對誤差小于文獻(xiàn)[9]和文獻(xiàn)[10]；當(dāng)網(wǎng)格比λ=1/2時，當(dāng)網(wǎng)格數(shù)N=16時，文獻(xiàn)[9]的計(jì)算結(jié)果發(fā)散；本文格式的最大絕對誤差要比文獻(xiàn)[10]小5～8個數(shù)量級，而且文獻(xiàn)[10]的方法只具有二階精度，但本文格式仍具有四階精度；表4給出了網(wǎng)格比λ=3/8和λ=1/2時，N=10,n=200(其中n=t/τ為總的時間步數(shù))時的問題2的數(shù)值解和精確解，從表中可以看出本文格式的計(jì)算結(jié)果較文獻(xiàn)[9]和文獻(xiàn)[10]的計(jì)算結(jié)果更加精確，這有效地驗(yàn)證了本文格式的精確性和可靠性.當(dāng)λ=1/2時，從表4中的計(jì)算結(jié)果可以看出，文獻(xiàn)[9]的計(jì)算結(jié)果是收斂的，這說明文獻(xiàn)[9]的方法對所求解的問題有依賴性，取同樣的參數(shù)，計(jì)算問題2時結(jié)果收斂，而計(jì)算問題1時結(jié)果卻發(fā)散.

表3 　問題2在t=1時刻的最大絕對誤差及收斂階

表4　問題2在N=10的數(shù)值解和精確解

問題3是非齊次的擴(kuò)散方程，由于文獻(xiàn)[9]和文獻(xiàn)[10]中未給出關(guān)于非齊次問題的求解格式，因此，無法與文獻(xiàn)[9]和文獻(xiàn)[10]作比較.對于這個問題，與文獻(xiàn)[2]的高精度隱格式進(jìn)行了比較.問題3的邊界條件是關(guān)于時間t的函數(shù)，從表5中的計(jì)算結(jié)果可知，本文格式對此類問題是有效的，特別是對邊界的處理方法是可靠的；表5給出了問題3在λ=1/3和λ=1/2時的最大絕對誤差、CPU時間T及收斂階，由最大計(jì)算誤差可知，本文格式比文獻(xiàn)[2]更精確，從CPU時間T來看，本文格式計(jì)算所需的CPU時間T比文獻(xiàn)[2]的CPU時間少，特別是當(dāng)λ=1/2時，取時間t=5，網(wǎng)格數(shù)為N=256時，文獻(xiàn)[2]所用CPU時間T為117.0，而本文格式僅為63.15，這說明本文顯格式在計(jì)算網(wǎng)格數(shù)比較大且時間歩循環(huán)比較多的問題時具有一定的優(yōu)勢，在網(wǎng)格充分細(xì)化的計(jì)算過程中，能有效地節(jié)省存儲空間和降低計(jì)算成本.

表5　問題3在t=1時刻的最大絕對誤差、CPU時間及收斂階Rate

4　結(jié)論

本文對一維擴(kuò)散方程初邊值問題的顯式高精度緊致差分格式進(jìn)行了研究，得出以下結(jié)論：

1)采用四階緊致差分公式對空間二階導(dǎo)數(shù)進(jìn)行離散，采用泰勒級數(shù)展開的方法對時間導(dǎo)數(shù)項(xiàng)進(jìn)行離散，推導(dǎo)得到了一種數(shù)值求解一維擴(kuò)散方程的顯式高精度緊致差分格式，截?cái)嗾`差為O(τ2+τh2+h4),并通過Fourier方法分析了格式的穩(wěn)定性，其穩(wěn)定性條件為λ≤1/2.通過數(shù)值實(shí)驗(yàn)，驗(yàn)證了其精確性和穩(wěn)定性.

2)目前已經(jīng)發(fā)展的求解擴(kuò)散問題的高精度格式多為3層格式，如文獻(xiàn)[9-10,12-13]，對3層格式的計(jì)算，需要采用另外的格式進(jìn)行啟動，即完成第1個時間步的計(jì)算.而本文格式為2層格式，直接由初值步開始依次計(jì)算，并且本文格式與文獻(xiàn)[9-10]中的格式具有相同的穩(wěn)定性條件，但是本文格式的計(jì)算結(jié)果卻更為精確.

3)將本文格式與同階精度的隱式方法[2]進(jìn)行了比較，盡管本文方法的穩(wěn)定性不及隱式方法，但是在本文方法的穩(wěn)定性條件滿足的情況下，取相同的計(jì)算參數(shù)，本文方法可以得到較隱式方法[2]更為精確的計(jì)算結(jié)果，并且由于本文方法是顯式計(jì)算，因此較隱式方法的計(jì)算量大為減少，提高了問題的求解效率，節(jié)約了CPU時間.

[1]田振夫.一維對流擴(kuò)散方程的四階精度有限差分法[J].寧夏大學(xué)學(xué)報(bào)(自然科學(xué)版),1995,16(1):30-35.

TIAN Zhenfu.Fourth-order finite difference methods for one-dimensional convection-diffusion equations[J].Journal of Ningxia University(Natural Science Edition)，1995,16(1):30-35.

[2]葛永斌,田振夫,詹詠，等.求解擴(kuò)散方程的一種高精度隱式差分方法[J].上海理工大學(xué)學(xué)報(bào),2005,27(2):107-112.DOI:10.3969/j.issn.1007-6735.2005.02.004.

GE Yongbin,TIAN Zhenfu,ZHAN Yong,et al.High order implicit difference method for the diffusion equation[J].Journal of University of Shanghai for Science and Technology,2005,27(2):107-112.DOI:10.3969/j.issn.1007-6735.2005.02.004.

[3]詹涌強(qiáng),張傳林.解拋物型方程的一族高精度隱式差分格式[J].應(yīng)用數(shù)學(xué)和力學(xué),2014,35(7):790-797.DOI:10.3879/j.issn.1000-0887.2014.07.008.

ZHAN Yongqiang,ZHANG Chuanlin.A family of high accuracy implicit difference schemes for solving parabolic equatons[J].Applied Mathematics and Mechanics,2014,35(7):790-797.DOI:10.3879/j.issn.1000-0887.2014.07.008.

[4]詹涌強(qiáng),張傳林.解拋物型方程的一個新的高精度隱格式[J].華中師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,48(2):168-170.

ZHAN Yongqiang,ZHANG Chuanlin.A new high accuracy implicit difference schemes for solving parabolic equatons[J].Journal of Huazhong Normal University(Natural Science Edition),2014,48(2):168-170.

[5]周敏,高學(xué)軍,董超.解拋物型方程的八點(diǎn)隱式差分格式[J].廣東工業(yè)大學(xué)學(xué)報(bào),2014,31(4):71-78.DOI:10.3969/j.issn.1007-7162.2014.04.013.

ZHOU Min,GAO Xuejun,DONG Chao.The implicitght difference scheme of eight points for solving the parabolic equtions[J].Journal of Guangdong University of Technology,2014,31(4):71-78.DOI:10.3969/j.issn.1007-7162.2014.04.013.

[6]于德浩,湯華中.偏微分方程數(shù)值解法[M].北京:科學(xué)出版社,2004:106-109.

YU Dehao,TANG Huazhong.Numerical solution of partial differential equation[M].Beijing:Science Press,2004:106-109.

[7]楊情民.解拋物型方程的一族顯格式[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào).1981,4:306-317.

YANG Qingmin.A family of explicit schemes for solving the parabolic equations[J].Journal of High School Mathematics,1981,4:306-317.

[8]曾文平.多維拋物型方程的分支絕對穩(wěn)定的顯式格式[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),1997,19(2):112-121.

ZENG Wenping.Explicit scheme for the absolute stability of the multi dimensional parabolic equations[J].Journal of High School Mathematics,1997,19(2):112-121.

[9]金承日.解拋物型方程的高精度顯式差分格式[J].計(jì)算數(shù)學(xué),1991,13(1):38-44.DOI:10.3321/j.issn:0254-7791.1991.01.001.

JIN Chengri.High accuracy explicit difference scheme for solving parabolic equations[J].Computational Mathematics,1991,13(1):38-44.DOI:10.3321/j.issn:0254-7791.1991.01.001.

[10]馬明書.一維拋物型方程的一個新的高精度顯式差分格式[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2001,22(2):156-160.DOI:10.3969/j.issn.1000-3266.2001.02.010.

MA Mingshu.A new high accuracy explicit difference scheme with branching stable for solving parabolic equation of one-dimension[J].Numerical Computation and Computer Application,2001,22(2):156-160.DOI:10.3969/j.issn.1000-3266.2001.02.010.

[11]馬明書.拋物型方程的一個新的顯格式[J].紡織高?；A(chǔ)科學(xué)學(xué)報(bào),2001,14(2):133-135.DOI:10.3969/j.issn.1006-8341.2001.02.010.

MA Mingshu.A new explicit difference scheme with high accuracy for solving parabolic type equation[J].Journal of Basic Science of Textile University,2001,14(2):133-135.DOI:10.3969/j.issn.1006-8341.2001.02.010.

[12]徐金平,單雙榮.解拋物型方程的一個高精度顯式差分格式[J].華僑大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,30(4):473-475.DOI:10.11830/ISSN.1000-5013.2009.04.0473.

XU Jinping,SHAN Shuangrong.An explicit difference scheme with high-order accuracy for solving parabolic equation[J].Journal of Huaqiao University(Natural Science),2009,30(4):473-475.DOI:10.11830/ISSN.1000-5013.2009.04.0473.

[13]袁權(quán)龍,詹再東.拋物方程高精度高穩(wěn)定顯格式研究[J].山西師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,23(4):27-30.DOI:10.3969/j.issn.1009-4490.2009.04.007.

YUAN Quanlong,ZHAN Zaidong.High accuracy and stable explicit scheme for the parabolic equation[J].Journal of Shanxi Normal University(Natural Science Edition),2009,23(4):27-30.DOI:10.3969/j.issn.1009-4490.2009.04.007.

[14]WEN ruihao,Shao Hongzhu.Domain decomposition schemes with high-order accuracy and unconditional stability[J].Appl Math Comput,2013,219:6170-6181.DOI:http://dx.doi.org/10.1016/j.amc.2012.12.001.

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

(責(zé)任編輯：王蘭英)

High accuracy explicit compact difference scheme for the diffusion equation

YANG Xiaojia,WANG Yan

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

Based on the fourth-order compact difference formula of the second-order derivative in spatial direction and Taylor series expansion in time direction,a high-order explicit compact difference scheme for the one dimentional diffusion equation is developed.The stability of the scheme is analyzed by Fourier method,and the condition of stability isλ≤1/2 (λis the mesh ratio).Numerical experiments are carried out to verify the accuracy and reliability of the present scheme.

diffusion equation;high accuracy;compact scheme;explicit scheme;finite difference method

10.3969/j.issn.1000-1565.2016.02.002

2015-09-08

寧夏高等學(xué)?？茖W(xué)技術(shù)研究項(xiàng)目(NGY2013019)

楊曉佳(1988—)，女，寧夏吳忠人，寧夏大學(xué)在讀碩士研究生.E-mail：yang_xiaoj@sina.com

王燕(1976—)，女，寧夏銀川人，寧夏大學(xué)副教授，主要從事偏微分方程數(shù)值解法的研究.

E-mail：wangy@nxu.edu.cn

O241.82

A

1000-1565(2016)02-0117-07