求解波動(dòng)方程的2種顯式高精度緊致差分格式

姜蘊(yùn)芝,葛永斌

(寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院,寧夏 銀川 750021)

求解波動(dòng)方程的2種顯式高精度緊致差分格式

姜蘊(yùn)芝,葛永斌*

(寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院,寧夏 銀川 750021)

針對(duì)一維波動(dòng)方程，空間采用四階Padé逼近，時(shí)間采用中心差分離散得到了一種時(shí)間二階、空間四階精度的顯式緊致差分格式，其截?cái)嗾`差為O(τ2+h4).之后采用截?cái)嗾`差余項(xiàng)修正的方法對(duì)時(shí)間離散進(jìn)行改進(jìn)，改進(jìn)后的格式的截?cái)嗾`差為O(τ4+τ2h2+h4)，即格式具有整體四階精度.然后，通過Fourier方法分析了2種格式的穩(wěn)定性.最后，通過數(shù)值實(shí)驗(yàn)驗(yàn)證了本格式的精確性和可靠性.

波動(dòng)方程； Padé逼近； 緊致格式； 顯式差分； 穩(wěn)定性

波動(dòng)方程是一類重要的雙曲型偏微分方程，在數(shù)學(xué)、物理、化學(xué)等領(lǐng)域內(nèi)有著廣泛的應(yīng)用[1].對(duì)這類方程進(jìn)行數(shù)值求解的方法主要有有限差分法、有限元法和有限體積法等.就有限差分法而言，對(duì)該類問題研究的理論成果有：文獻(xiàn)[2]采用二階中心差分格式和非均勻網(wǎng)格離散，提出了一種求解一維波動(dòng)方程在非均勻時(shí)間網(wǎng)格上的三層顯式差分格式，該格式具有一階精度；文獻(xiàn)[3]利用泰勒級(jí)數(shù)展開及待定系數(shù)法建立了一種求解一維波動(dòng)方程的三層隱格式，該格式是條件穩(wěn)定的，并且具有四階精度；文獻(xiàn)[4]將Runge-Kutta方法應(yīng)用到哈密頓系統(tǒng)中并與辛格式相結(jié)合，提出了求解一維波動(dòng)方程的一類顯式辛方法，該方法具有二階精度，并且是條件穩(wěn)定的；文獻(xiàn)[5]采用三次樣條公式推導(dǎo)出精度分別為O(τ2+h2)、O(τ2+h4)、O(τ4+h2)和O(τ4+h4)的4種差分格式，并采用Fourier方法分析了格式的穩(wěn)定性；文獻(xiàn)[6]對(duì)一維二階波動(dòng)方程提出了具有二階精度的精細(xì)時(shí)程積分法，該方法能夠在時(shí)間方向上精確計(jì)算，空間方向的局部截?cái)嗾`差為O(h2)，并且該方法是無條件穩(wěn)定的；文獻(xiàn)[7]通過加權(quán)平均和緊致差分離散的思想提出了2種精度分別為O(τ2+h4)和O(τ4+h4)的隱式緊致差分格式；文獻(xiàn)[8]在Crank-Nicolson格式的基礎(chǔ)上設(shè)計(jì)了重疊型區(qū)域分解的并行算法，該算法的最優(yōu)逼近階為O(τ2+h2)；文獻(xiàn)[9]利用三次樣條插值，提出了求解一維波動(dòng)方程的一種三層隱式差分格式，該格式最優(yōu)能夠達(dá)到時(shí)間二階，空間四階精度，并且是無條件穩(wěn)定的；文獻(xiàn)[10]通過四次樣條函數(shù)與廣義梯形算法相結(jié)合的方法提出了一維波動(dòng)方程的一類兩層差分格式，其精度為O(τ2+h4)，當(dāng)選擇適當(dāng)?shù)膮?shù)時(shí)，其精度可提高到O(τ3+h4)；文獻(xiàn)[11]采用四階緊致差商逼近公式及加權(quán)平均思想，提出了2種精度分別為O(τ2+h4)和O(τ4+h4)的交替方向隱式格式，前者是無條件穩(wěn)定的，后者是條件穩(wěn)定的；文獻(xiàn)[12]提出了一種求解一維波動(dòng)方程的高精度隱式差分格式，該格式是無條件穩(wěn)定的，并且具有時(shí)間二階、空間四階精度.本文將建立2種顯式緊致差分格式，為此，考慮如下一維波動(dòng)方程的初邊值問題:

(1)

(2)

(3)

其中,u(x,t)是待求未知函數(shù)，a為波動(dòng)系數(shù)，f(x,t)為非齊次項(xiàng)，φ(x)、ψ(x)、g0(t)、g1(t)為已知函數(shù)，且具有充分的光滑性.

1 顯式緊致差分格式的構(gòu)造

1.1 CTFS格式對(duì)初始時(shí)間層的離散.利用泰勒展開公式有

(4)

將(2)式代入(4)式中，進(jìn)行整理并舍去其截?cái)嗾`差項(xiàng)有

(5)

(6)

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

(7)

(8)

(9)

對(duì)上式進(jìn)行整理，并舍去其截?cái)嗾`差項(xiàng)有

(10)

1.2 FTFS格式為了使時(shí)間精度與空間精度能夠相匹配，使格式整體精度達(dá)到四階，下面對(duì)上述差分格式進(jìn)行改進(jìn)，進(jìn)一步提高時(shí)間精度，為此，對(duì)初始時(shí)間層的離散采用文獻(xiàn)[7]中的方法，有

(11)

其中λ=τ/h.對(duì)于時(shí)間二階導(dǎo)數(shù)項(xiàng)的離散，采用中心差分格式離散后保留其截?cái)嗾`差主項(xiàng)，可得

(12)

又由(1)式可得

(13)

將上式代入(12)式中有

(14)

對(duì)(14)式進(jìn)行整理，并舍去其截?cái)嗾`差項(xiàng)有

又由(1)式，對(duì)上式進(jìn)行變形有

(15)

(15)式即為整體四階精度的顯式緊致差分格式，記為FTFS格式.由格式的構(gòu)造過程可知，該差分格式的截?cái)嗾`差為O(τ4+τ2h2+h4).

2 穩(wěn)定性分析

引理 1[1]實(shí)系數(shù)二次方程μ2-bμ-c=0的根按模不大于1的充分必要條件為|b|≤1-c≤2.

對(duì)于(6)式有

(16)

對(duì)上式進(jìn)行化簡(jiǎn)整理有

(17)

(18)

令Uj=(uj,vj)T，并將(17)式代入進(jìn)行整理有

(19)

從而可得CTFS格式(10)的誤差增長(zhǎng)矩陣為

(20)

令λ=τ/h，r=aλ,則得上述誤差增長(zhǎng)矩陣的特征方程為

(21)

(22)

(23)

(24)

將(17)式代入(24)式進(jìn)行整理，可得FTFS格式的誤差增長(zhǎng)矩陣為

(25)

則誤差增長(zhǎng)矩陣的特征方程為

(26)

其中

因此，可得該格式穩(wěn)定的充要條件為

(27)

上式等價(jià)于如下2個(gè)不等式:

(28)

對(duì)于不等式(28)易得

(30)

由于不等式對(duì)任意ω的取值都要成立，所以有

(31)

對(duì)于不等式(29)，由求根公式易得其等式解為

(32)

(33)

3 數(shù)值算例

為了驗(yàn)證本文所提2種格式的精確性和可靠性，現(xiàn)考慮如下2個(gè)具有精確解的初邊值問題.

問題 1[7]:

其精確解為u(x,t)=sin(πt)sin(πx).

問題 2 :

其精確解為u(x,t)=te-πtsin(πx).

表1～3給出了問題1的計(jì)算結(jié)果.表1采用本文CTFS格式與文獻(xiàn)[7]中的四階格式進(jìn)行了計(jì)算.由于這2種格式的精度均為時(shí)間二階、空間四階，因此取τ=2h2，計(jì)算了t=0.5時(shí)刻取不同h時(shí)(τ也相應(yīng)不同)的L∞和L2范數(shù)誤差和收斂階.eL∞和eL2范數(shù)及收斂階的定義如下：

(34)

(35)

(36)

由表1的數(shù)據(jù)可以看出，2種格式空間均達(dá)到了四階精度，而本文的CTFS格式要比文獻(xiàn)[7]中的四階格式更為精確.表2給出了本文FTFS格式和文獻(xiàn)[7]中的四階格式的計(jì)算結(jié)果，由于2種格式均具有整體四階精度，因此取τ=h，計(jì)算了t=0.5時(shí)刻取不同h時(shí)的L∞和L2范數(shù)誤差和收斂階.可以看出，2種格式幾乎具有相同的精度.

因此，當(dāng)|a|λ大過1之后，計(jì)算結(jié)果是發(fā)散的.

表1 問題1當(dāng)τ=2h2時(shí)，在t=0.5時(shí)刻對(duì)不同h的L2和L∞范數(shù)誤差及收斂階Table 1 The L2,L∞ norm error and the convergence rate with τ=2h2 at t=0.5 for different h for Problem 1

注:1.46E-3=1.46×10-3，下同.

表2 問題1當(dāng)τ=h時(shí)，在t=0.5時(shí)刻對(duì)不同h的L2和L∞范數(shù)誤差及收斂階Table 2 The L2,L∞ norm error and the convergence rate with τ=h at t=0.5 for different h for Problem 1

表3 問題1當(dāng)h=1/32時(shí)，對(duì)不同λ的L∞范數(shù)誤差Table 3 The L∞ norm error with h=1/32 for different λ for Problem 1

表4～6給出了問題2的計(jì)算結(jié)果.

表4 問題2當(dāng)τ=h2時(shí)，在t=1時(shí)刻對(duì)不同h的L2和L∞范數(shù)誤差及收斂階Table 4 The L2,L∞ norm error and the convergence rate with τ=h2 at t=1 for different h for Problem 2

表5 問題2當(dāng)τ=h時(shí)，在t=1時(shí)刻對(duì)不同h的L2和L∞范數(shù)誤差及收斂階Table 5 The L2,L∞ norm error and the convergence rate with τ=h at t=1 for different h for Problem 2

表6 問題2當(dāng)h=1/32時(shí)，對(duì)不同λ的L∞范數(shù)誤差Table 6 The L∞ norm error with h=1/32 for different λ for Problem 2

4 結(jié)論

致謝 寧夏大學(xué)自然科學(xué)基金項(xiàng)目(ZR1407)和寧夏大學(xué)研究生創(chuàng)新項(xiàng)目(GIP2016032)對(duì)本文給予了資助,謹(jǐn)致謝意.

[1] 戴嘉尊,邱建賢.偏微分方程數(shù)值解法[M].南京:東南大學(xué)出版社,2002.

[2] CHO C H.Stability for the finite difference schemes of the linear wave equation with nonuniform time meshes[J].Num Meth Part Diff Eqns,2013,29(3):1031-1042.

[3] 張?zhí)斓?孫傳灼.關(guān)于波動(dòng)方程的差分格式[J].山東工業(yè)大學(xué)學(xué)報(bào),1995,25(3):283-287.

[4] 孫耿.波動(dòng)方程的一類顯式辛格式[J].計(jì)算數(shù)學(xué),1997,19(1):1-10.

[5] 齊遠(yuǎn)節(jié),劉利斌.求解二階波動(dòng)方程的三次樣條差分方法[J].大學(xué)數(shù)學(xué),2011,27(1):59-64.

[6] 金承日,呂萬金.二階雙曲型方程的精細(xì)時(shí)程積分法[J].計(jì)算力學(xué)學(xué)報(bào),2003,20(1):113-115.

[7] 葛永斌,朱琳,田振夫.求解波動(dòng)方程的高精度緊致隱式差分方法[J].寧夏大學(xué)學(xué)報(bào)(自然科學(xué)版),2005,26(4):297-299.

[8] 田敏,羊丹平.波動(dòng)方程的重疊型區(qū)域分解并行有限差分算法[J].山東大學(xué)學(xué)報(bào)(理學(xué)版),2007,42(2):28-38.

[9] RASHIDINIA J,JALILIAN R,KAZEMI V.Spline methods for the solutions of hyperbolic equations[J].Appl Math Comput,2007,190(1):882-886.

[10] LIU T,LIU L,XU H,et al.A new two level difference scheme for solving one-dimensional second-order hyperbolic equations[C]//2012 Fifth International Joint Conference on Computational Sciences and Optimization,2010:218-221.

[11] 馬月珍,李小剛,葛永斌.二維波動(dòng)方程的高精度交替方向隱式方法[J].四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,33(2):179-183.

[12] 崔進(jìn),吳宏偉.一類波動(dòng)方程初邊值問題的高階差分格式[J].應(yīng)用數(shù)學(xué),2014,27(1):166-174.

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

[14] 葛永斌,吳文權(quán),田振夫.二維波動(dòng)方程的高精度隱格式及其多重網(wǎng)格算法[J].廈門大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,42(6):691-696.

[15] 余躍玉.一種Caputo型時(shí)間分?jǐn)?shù)階波動(dòng)方程的差分方法[J].四川師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,37(4):524-528.

2010 MSC：35A35; 65M99

(編輯 陶志寧)

Two Kinds of Explicit High Order Compact Difference Schemes for Solving Wave Equations

JIANG Yunzhi,GE Yongbin

(CollegeofSchoolofMathematicalandStatistics,NingxiaUniversity,Yinchuan750021,Ningxia)

In this paper,an explicit compact difference scheme is obtained for solving the one dimensional wave equation.The truncation error of the scheme isO(τ2+h4).It’s constructed by applying the fourth-order accurate Padé approximation in space and the second-order accurate central difference in time.Then,the remainder of the truncation error correction method is employed to improve the accuracy of the discretization of time,the truncation error of the improved scheme isO(τ4+τ2h2+h4),which means the scheme has an overall fourth-order accuracy.And then,the stability conditions of the two schemes are obtained by the Fourier method.Finally,the accuracy and the reliability of the present two schemes are verified by numerical experiments.

wave equation; Padéapproximation; compact scheme; explicit difference; stability

2015-01-16

國(guó)家自然科學(xué)基金(11361045)

O242.1

A

1001-8395(2017)02-0177-07

10.3969/j.issn.1001-8395.2017.02.006

*通信作者簡(jiǎn)介：葛永斌(1975—),男,教授,主要從事偏微分方程數(shù)值解法和計(jì)算流體力學(xué)的研究,E-mail:gyb@nxu.edu.cn