二維波動方程的高階精度緊致顯式差分格式及穩(wěn)定性分析

摘 要：針對具有初邊值問題的二維波動方程，提出了一種數(shù)值求解該方程的高階緊致顯式有限差分格式。首先，根據(jù)相關(guān)文獻對導數(shù)的離散近似，得到周期邊界條件下的六階緊致差分格式。其次，在空間方向上，邊界節(jié)點導數(shù)項利用原方程代入的方法進行計算，而內(nèi)部節(jié)點的導數(shù)項利用六階緊致差分公式近似，使空間精度達到六階。同時，在時間方向上，利用泰勒級數(shù)展開公式、原方程代入以及中心差分公式推導出時間層的二階精度差分格式，為了將整體上的時間精度由二階提高至四階，采用外推算法實現(xiàn)時間層的高階近似。再次，再利用傅里葉分析法對該格式的穩(wěn)定性進行分析，得到在此精度下的穩(wěn)定性條件，即|a|λ∈［0，1276］。最后，通過數(shù)值實驗驗證了所提出的HOCE（6，4）格式的高效性和準確性。

關(guān)鍵詞：波動方程；中心差分；緊致差分；外推算法；穩(wěn)定性

DOI：10.15938/j.jhust.2024.03.017

中圖分類號： O241

文獻標志碼： A

文章編號： 1007-2683（2024）03-0141-08

High Order Compact Explicit Difference Schemes and

Stability Analysis for Two-dimensional Wave Equations

SUN Yang1， SONG Linlin1， AI Xiaohui2

（1.School of Science， Harbin University of Science and Technology， Harbin 150080， China；

2.School of Science， Northeast Forestry University， Harbin 150040， China）

Abstract：In this paper， a high-order compact explicit finite difference scheme is proposed to numerically solve two-dimensional wave equations with initial boundary value problems. First of all， according to the discrete approximation of derivatives in the existing literature， a sixth-order compact difference scheme with periodic boundary conditions is obtained. Then， in the spatial direction， the derivative term of the boundary node is calculated by substituting the original equation， and the derivative term of the internal node is approximated by the sixth order compact difference formula， so that the spatial accuracy can reach the sixth order. For the time direction， the second-order accuracy difference scheme of the time layer is derived by using the Taylor series expansion formula， the original equation substitution and the central difference formula. In order to improve the overall time accuracy from the second order to the fourth order， the Richardson extrapolation method is used to realize the high-order approximation of the time layer. Furthermore the stability of the scheme is analyzed by Fourier analysis method and the stability condition is |a|λ∈［0，1276］. Finally， the efficiency and accuracy of the proposed HOCE（6，4） scheme are verified by numerical experiments.

Keywords：wave equation; central difference; compact finite difference; Richardson extrapolation; stability

0 引 言

波動方程作為一類重要的描述波動現(xiàn)象的雙曲型偏微分方程，一直以來受到人們的廣泛關(guān)注。其在彈性力學、流體力學、電磁學及氣象學等諸多自然科學領(lǐng)域都有著較為深入的研究與應(yīng)用［1］。眾所周知，實際問題具有高度的復雜性，求解偏微分方程的解析解有時較為困難。因此，對于偏微分方程數(shù)值解法的研究具有重要意義。

目前，數(shù)值求解方法主要有：有限差分法、有限元法、有限體積法以及譜方法等。其中以離散點的差商近似導數(shù)作為研究手段的有限差分法具有形式簡單、易于編程實現(xiàn)和計算精度高等優(yōu)點，是求解偏微分方程的有效方法［2］。有限差分法的產(chǎn)生時間較早，其數(shù)學理論也已經(jīng)較為成熟，常見的分析差分法穩(wěn)定性的方法主要有Von Neumann方法、矩陣分析方法和能量分析方法［3-4］。

低精度的差分離散格式雖然計算相對穩(wěn)定，但是由于其較低的計算精度和分辨率，往往導致數(shù)值計算失真，因此眾多的計算數(shù)學工作者專注于高精度、高分辨率的計算格式研究。研究表明，高精度格式較低精度格式可以有效地減小對步長的限制，即使所使用的網(wǎng)格為粗網(wǎng)格，也能夠獲得低精度格式在較細網(wǎng)格上的計算精度，即在計算域內(nèi)利用比較少的網(wǎng)格點就可以獲得較高的計算精度，從而提高計算效率。此外，有研究表明，高精度的數(shù)值方法可以有效抑制數(shù)值色散。然而，隨著計算精度的提高，高精度格式可能會出現(xiàn)數(shù)值不穩(wěn)定。因此，發(fā)展波動方程的高精度（高于二階精度）且穩(wěn)定性好的差分格式是重要且具有實際意義的研究課題［5-6］。

近年來，較多學者對波動方程進行了大量的研究，并提出了很多精確有效的方法。其中，文［7］針對一維波動方程，應(yīng)用Runge-Kutta方法，提出了一種具有嚴格的穩(wěn)定性條件的的二階精度的顯式格式。在文［8-10］中利用Richardson外推法，將時間二階精度提高至四階，但需要計算細網(wǎng)格的二階精度解進一步得到粗網(wǎng)格上的四階精度解，從而降低計算效率。綜上，許多學者在解決此類問題時，使用的方法在時間或空間上精度較低［11-12］。對于多維問題，高精度差分格式因其較小的數(shù)值耗散等優(yōu)點有著廣泛的應(yīng)用，但是也存在明顯的問題：即在處理邊界條件問題上存在不足，所需精度越高節(jié)點就越多，導致處理邊界時計算區(qū)域兩端會出現(xiàn)所謂“冒點”的情況，從而增加邊界處理難度，同時破壞了系數(shù)矩陣的對稱性、增加了計算時間，導致高精度格式的計算效率有所下降。如果采用隱式方法，由于有限差分算子的有理函數(shù)將替代差分算子，進而每個時間步都要解決三對角系統(tǒng)，使計算量增加。常用的解決方法是交替方向隱（alternating direction implicit， ADI）格式和局部一維（local one-dimension， LOD）格式［13-14］。ADI格式是將高維問題轉(zhuǎn)化為一維問題進行計算，從而減少計算量。在引入ADI格式用于求解雙曲型方程后，有很多學者將其推廣到二維非線性雙曲方程，在文［15］中有詳細介紹。但是，ADI格式對計算區(qū)域有限制條件。LOD格式相比ADI格式解決了這一問題。在文［16］中，討論了基于Richardson外推的時間和空間四階LOD格式。在文［17］中提出包含LOD的三層隱式格式。然而，以上提出的這些方法都是非緊致的，或者只有在特殊條件下才能使用。為了解決這些問題，一些學者致力于各種緊致的高階差分格式的研究與應(yīng)用。高階緊致差分格式是利用不同節(jié)點處空間導數(shù)值的線性組合等于不同節(jié)點處函數(shù)的線性組合進行差分，然后根據(jù)不同的精度要求，通過待定系數(shù)法得到這些組合系數(shù)。近年來，相關(guān)數(shù)值實驗可以證明，相比非緊致和低階方法，高階緊致差分格式要求的計算節(jié)點相對較少、邊界處理相對簡單且在一定程度上具有更好的分辨率，并且能夠證明緊致格式的截斷誤差通常比非緊致格式的截斷誤差在同階格式下小四到六倍［18-19］。

本文針對二維波動方程，應(yīng)用有限差分思想，構(gòu)造一種高階精度緊致顯C2g2dZTCI14Ktlqq0X8bPA==式有限差分格式。由于本文格式是顯式的，因此不需要任何迭代過程，能夠有效的縮短計算時間。以往數(shù)值結(jié)果表明，緊致格式優(yōu)于相應(yīng)的同階顯式格式。因此，本文格式是將常見的二階或四階空間精度提高到六階，時間方向上利用Richardson外推，將時間精度由二階提高至四階的高階緊致顯式差分格式。然后利用Fourier分析法分析該格式的穩(wěn)定性。為了方便起見，將本文的空間六階、時間四階計算格式簡寫為HOCE（6，4）格式。為驗證其計算效率和高精度特性，本文在數(shù)值實驗中選擇了其它兩種計算格式：空間四階、時間二階（簡稱：HOCE（4，2））格式［21］和空間六階、時間二階（簡稱：HOCE（6，2））格式，并且與本文HOCE（6，4）格式的計算結(jié)果進行比較分析，從而驗證了本文格式的高精度特性以及在計算效率方面的較強優(yōu)勢。

1 問題描述

對于如下二維波動方程的初邊值問題

2ut2=a2（2ux2+2uy2）+f（x，y，t）

（x，y，t）∈D×（0，T］

u（x，y，0）=φ（x，y）

u（x，y，0）t=ψ（x，y）

u（0，y，t）=g0（y，t），u（1，y，t）=g1（y，t）

u（x，0，t）=h0（x，t），u（x，1，t）=h1（x，t）（1）

由于進行數(shù)值求解，需要增加的周期邊界條件如下：

u（x+1，y，t）=－u（x，y，t）

u（x，y+1，t）=－u（x，y，t）（2）

式中：u（x，y，t）是未知函數(shù)；D=［0，1］×［0，1］；f（x，y，t）、φ（x，y）、ψ（x，y）、g0（y，t）、g1（y，t）、h0（x，t）、h1（x，t）均為已知函數(shù)，且充分光滑；a為波動系數(shù)。

在計算區(qū)域內(nèi)，進行等距網(wǎng)格剖分，h是空間步長，τ是時間步長，（xi，yj，tn）是網(wǎng)格節(jié)點，其中，xi=ih，yj=jh，tn=nτ，i，j=0，1，…，N，n=0，1，…，M，用uni，j近似表示（xi，yj，tn）的近似值。

2 HOCE（6，4）格式

首先，根據(jù)文［18-19］中的二階導數(shù)的離散，可以得到如下公式

αu″j－1+u″j+αu″j+1=

a1uj+1－2uj+uj－1h2+a2uj+2－2uj+uj－24h2（3）

式中α，a1，a2為自由參數(shù)，且決定格式的精度，在這里可取

a1=43（1－α），a2=13（10α－1）。

當α=211時，即a1=1211， a2=311，此時式（3）即為六階精度。在周期邊界條件下，將上述參數(shù)α，a1，a2代入式（3），可得如下的六階精度格式

211u″j－1+u″j+211u″j+1=

1211uj+1－2uj+uj－1h2+311uj+2－2uj+uj－24h2

即

8u″j－1+44u″j+8u″j+1=

3uj+2+48uj+1－102uj+3uj－2+48uj－1h2（4）

當參數(shù)α=110，a1=1210，a2=0時，式（3）為四階精度，即

u″j－1+10u″j+u″j+1=12uj+1－2uj+uj－1h2（5）

2.1 計算（2ux2）ni，j和（2uy2）ni，j的值

空間內(nèi)部節(jié)點利用六階緊致差分格式，則有

8（2ux2）ni+1，j+44（2ux2）ni，j+8（2ux2）ni-1，j=

3uni+2，j+48uni+1，j－102uni，j+3uni－2，j+48uni－1，jh2（6）

式中i=2，3，…，N－2，j=0，1，…，N。

8（2uy2ni，j+1+44（2uy2）ni，j+8（2uy2ni，j－1=

3uni，j+2+48uni，j+1－102uni，j+3uni，j－2+48uni，j－1h2（7）

式中：j=2，3，…，N－2；i=0，1，…，N。

對空間邊界點進行處理，由式（1）和式（2）可得

（2ux2）n0，j=1a2（2ut2-f）n0，j-（2uy2）n0，j=

1a2（2g0t2-f）n0，j－（2g0y2n0，j（8）

式中j=0，1，…，N。

（2ux2nN，j=1a2（2ut2-f）nN，j－（2uy2nN，j=

1a2（2g1t2-f）nN，j－（2g1y2nN，j（9）

式中j=0，1，…，N。

（2uy2）n0，j=1a2（2ut2-f）ni，0-（2uy2）ni，0=

1a2（2h0t2-f）ni，0－（2h0y2ni，0（10）

式中i=0，1，…，N。

（2uy2）ni，N=1a2（2ut2-f）ni，N-（2uy2）ni，N=

1a2（2h1t2-f）ni，N－（2h1y2ni，N（11）

式中i=0，1，…，N。

2.2 第一個時間層的離散

利用Taylor公式展開得到

u1i，j=u0i，j+τ（ut）0i，j+τ22（2ut2）0i，j+O（τ3）=

u0i，j+τ（ut0i，j+τ22（2ut2）0i，j+O（τ3）（12）

利用式（1），對上式右端項進行整理可得

（2ut2）0i，j=a2［（2ux20i，j+（2uy2）0i，j］+f0i，j=

a2［（φxx）i，j+（φyy）i，j］+f0i，j（13）

將式（13）代入式（12），化簡并舍去截斷誤差項得

u1i，j=φ（xi，yj）+τψ（xi，yj）+

τ22［a2φxx（xi，yj）+a2φyy（xi，yj）+f（xi，yj，0）］（14）

2.3 其他時間層

利用中心差分格式近似時間的二階導數(shù)，并利用Taylor公式展開，可得

un+1i，j－2uni，j+un－1i，jτ2=a2［（2ux2）ni，j+（2uy2）ni，j］+

fni，j+o（τ2）（15）

對式（15）化簡整理，并舍去截斷誤差項得到

un+1i，j=2uni，j－un－1i，j+a2τ2［（2ux2）ni，j+（2uy2）ni，j］+fni，j（16）

式（16）即為空間六階精度，時間二階精度HOCE/kyzO3GOmU/AWuTjHnRhFQ==（6，2）的顯式緊致差分格式，該格式的截斷誤差為

O（h6+τ2h2+τ2）。

2.4 時間層的外推解

式（16）的格式在空間上為六階精度，但時間精度僅具有二階，考慮采用如下的Richardson外推法將時間精度提高到四階

uni，j（τ，h）=（4u2ni，j（τ/2，h）－uni，j（τ，h））/3（17）

其中：式（17）左端項uni，j（τ，h）為所求的四階精度解，右端項uni，j（τ，h）為格式（16）在時間步長為τ時的解，u2ni，j（τ/2，h）為時間步長取τ/2時的解。

3 穩(wěn)定性分析

利用Fourier方法分析格式穩(wěn)定性。這里，假設(shè)f精確無誤差，且令uni，j=ξnei（σ1xi+σ2yj），（uxx）ni，j=ηnei（σ1xi+σ2yj），（uyy）ni，j=γnei（σ1xi+σ2yj），式中ξ、η、γ為振幅，σ1、σ2為波數(shù)，i=－1為虛數(shù)單位。

引理1［20］ 實系數(shù)二次方程λ2－bλ－c=0的根按模不大于1的充要條件為|b|≤1－c≤2。

由已知條件

uni+1，j=ξnei（σ1xi+1+σ2yj）=

ξnei［σ1（xi+h）+σ2yj］=ξnei（σ1xi+σ2yj）eiσ1h

同理

uni－1，j=ξnei（σ1xi+σ2yj）e－iσ1h

uni+2，j=ξnei（σ1xi+σ2yj）e2iσ1h，uni－2，j=ξnei（σ1xi+σ2yj）e－2iσ1h

uni，j+1=ξnei（σ1xi+σ2yj）eiσ2h，uni，j－1=ξnei（σ1xi+σ2yj）e－iσ2h

uni，j+2=ξnei（σ1xi+σ2yj）eiσ2h，uni，j－2=ξnei（σ1xi+σ2yj）e－iσ2h

且可得到

（uxx）ni+1，j，（uxx）ni-1，j，（uxx）ni+2，j，（uxx）ni-2，j

（uyy）ni，j+1，（uyy）ni，j-1，（uyy）ni，j+2，（uyy）ni，j-2

則將上式代入式（6）化簡整理得

ηnei（σ1xi+σ2yj）［16cosσ1h+44］=

12h2ξnei（σ1xi+σ2yj）（cosσ1h－1）（cosσ1h+9）（18）

即

ηn=3h2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11ξn（19）

同理，式（7）化簡為

γnei（σ1xi+σ2yj）［16cosσ2h+44］=

12h2ξnei（σ1xi+σ2yj）（cosσ2h－1）（cosσ2h+9）（20）

即

γn=3h2（cosσ2h－1）（cosσ2h+9）4cosσ2h+11ξn（21）

下面進行穩(wěn)定性分析，令vn+1ij=unij，f≡0，式（16）的矩陣形式如下

un+1ijvn+1ij=2－110un+1ijvnij+

a2τ2000（uxx）nij+（uyy）nij（vxx）nij+（vyy）nij（22）

令Unij=（unij，vnij）T，Unij=ξnei（σ1h+σ2h），并代入式（22），化簡為

ξn+1=2－110ξn+a2τ2000（ηn+γn）（23）

將式（19）與式（21）代入式（23），記λ=τh，r=aλ，令

C=2+3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11（24）

可得該格式的誤差增長矩陣為

D=C－110

則對應(yīng)的特征方程是

μ2－2+3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11μ+1=0（24）

由引理1可知，在上述方程式（24）中，

b=2+3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11

c=－1，則該格式穩(wěn)定的充要條件為

2+3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11≤2（25）

式（25）等價為以下兩個不等式

3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11≤0（26）

－4≤3r2（cosσ1h－1）（cosσ1h+9）4cosσ1h+11+

（cosσ2h－1）（cosσ2h+9）4cosσ2h+11（27）

對于式（26），不等式恒成立。

對于式（27），化簡為r2≤43A，式中

A=（1－cosσ1h）（cosσ1h+9）4cosσ1h+11+

（1－cosσ2h）（cosσ2h+9）4cosσ2h+11

令α=1－cosσ1h，β=cosσ2h+9，則

A=α（10－α）15－4α+β（10－β）15－4β

且α∈［0，2］，β∈［0，2］，經(jīng)過特殊取值可發(fā)現(xiàn)，當α=β=2時，有

r2≤724

則該格式的穩(wěn)定性條件為

|a|λ∈［0，1276］

4 數(shù)值算例

為驗證本文格式的高精度和高分辨率特性，選擇比較HOCE（4，2）格式［21］，HOCE（6，2）格式與本文格式?？紤]以下方程的初邊值問題，計算了不同空間步長的不同時刻的L∞和L2范數(shù)誤差以及計算時間。其中L∞和L2范數(shù)誤差的定義如下：

L∞-error=maxi，j|uni，j－u（xi，yj，tn）|

L2-error=h2∑i，j［uni，j－u（xi，yj，tn）］2

算例1［21］

2ut2=2ux2+2uy2

u（0，y，t）=0，u（1，y，t）=0

u（x，0，t）=0，u（x，1，t）=0

u（x，y，0）=sin（πx）sin（πy）

u（x，y，0）t=0（28）

其精確解為u（x，y，t）=sin（πx）sin（πy）cos（2πt）。

通過比較HOCE（4，2）格式和HOCE（6，2）格式及本文格式的計算結(jié)果。即當h=1/10，1/20，…，1/320時，計算了τ=0.0002時在第100個時間步上，不同空間步長的范數(shù)誤差。其中，表1為L∞范數(shù)誤差，表2為L2范數(shù)誤差。由表1和表2可知，當空間步長不斷減小時，L2和L∞范數(shù)誤差也相應(yīng)減小，并且本文的計算結(jié)果優(yōu)于HOCE（4，2）和HOCE（6，2）的格式。表3為通過數(shù)值算例進一步說明本文所提出的HOCE（6，4）格式在計算時間方面具有高效性。對于不同的h，由于格式的計算復雜度不同，HOCE（6，4）格式計算時間較其它兩種格式有一定的增加。然而，若我們加密網(wǎng)格，在空間步長為h2時，由計算結(jié)果能夠看出HOCE（6，4）格式相比其它格式的計算時間無明顯差別，而計算誤差大為減少，從而充分體現(xiàn)了本文HOCE（6，4）格式的高精度的優(yōu)點。

算例2［21］

2ut2=2ux2+2uy2

u（0，y，t）=0，u（1，y，t）=0

u（x，0，t）=0，u（x，1，t）=0

u（x，y，0）=0

u（x，y，0）t=2πsin（πx）sin（πy）（29）

其精確解為

u（x，y，t）=sin（2πt）sin（πx）sin（πy）

表4和表5給出了本文HOCE（6，4）、本文HOCE（6，2）以及文［21］中的HOCE（4，2）格式的計算結(jié)果。計算了τ=0.0025，在t=1時，不同空間步長的范數(shù)誤差。即，L∞+范數(shù)誤差和L2范數(shù)誤差?？芍?，當空間步長不斷減小時，范數(shù)誤差也相應(yīng)減小，并且本文的計算結(jié)果相對較好。由表6可以得到，本文格式的計算相對簡單。

5 結(jié) 論

本文針對二維波動方程，利用六階緊致格式計算二階導數(shù)。由于緊致格式具有基架少，分辨率高等諸多優(yōu)點，因此特別適用于波動方程等偏微分方程的高精度計算。時間方向上利用Richardson外推，將時間二階提高到四階精度。利用Fourier分析法求出了該格式的穩(wěn)定區(qū)域為|a|λ∈［0，1276］。數(shù)值實驗方面，為驗證本文HOCE（6，4）格式的有效性和高分辨率特性，與HOCE（6，2）格式和HOCE（4，2）格式的計算結(jié)果進行比較分析。結(jié)果表明，不論計算精度還是計算時間方面，本文的HOCE（6，4）格式都具有很高的計算精度和計算效率。HOCE（6，4）格式在計算時間上較HOCE（6，2）格式并沒有明顯增加，但是在相同時間和空間步長的條件下，本文的HOCE（6，4）格式的計算誤差遠遠小于HOCE（6，2）格式和HOCE（4，2）格式，而且從計算結(jié)果可以看出：在粗網(wǎng)格下，本文的HOCE（6，4）格式的數(shù)值誤差要優(yōu)于空間步長減半的HOCE（4，2）格式的數(shù)值誤差。這充分驗證了高階精度能夠顯著提高差分格式的計算效率，因此可以說本文對計算格式的改進是成功的。

下一步作者考慮將HOCE（6，4）格式推廣到三維波動方程以及其它偏微分方程的數(shù)值計算當中，可以預見將會有更多較好的成果呈現(xiàn)。

參 考 文 獻：

［1］ 張?zhí)斓拢?孫傳灼. 關(guān)于波動方程的差分格式［J］. 山東工業(yè)大學學報， 1995（3）： 283.

ZHANG Tiande， SUN Chuanzhuo. On The Difference Scheme of Wave Equation［J］. Journal of Shandong University of Technology， 1995（3）： 283.

［2］ 王希， 張爽， 胡勁松. 求解廣義BBM-KdV方程的守恒型有限差分方法［J］. 哈爾濱理工大學學報， 2022， 27（4）： 147.

WANG Xi， ZHANG Shuang， HU Jinsong. Conservative Finite Difference Method for Solving Generalized BBM-KdV Equation. Journal of Harbin University of Science and Technology ， 2022， 27（4）： 147.

［3］ 冉茂華. 幾類分數(shù)階偏微分方程的有限差分方法［D］. 武漢. 華中科技大學， 2016.

［4］ 姜蘊芝， 葛永斌. 求解波動方程的2種顯式高精度緊致差分格式［J］. 四川師范大學學報：自然科學版， 2017， 40（2）： 177.

JIANG Yunzhi， GE Yongbin. Two Kinds of Explicit High Order Compact Difference Schemes for Solving Wave Equations［J］. Journal of Sichuan Normal University： Natural Science， 2017， 40（2）： 177.

［5］ HOU Baihui， DONG Liang. The Conservative Time High-Order AVF Compact Finite Difference Schemes for Two-Dimensional Variable Coefficient Acoustic Wave Equations［J］. Journal of Scientific Computing， 2019， 80（2）： 1279.

［6］ BRITT S， TURKEL E， TSYNKOV S. A High Order Compact Time/Space Finite Difference Scheme for The Wave Equation with Variable Speed of Sound［J］. Journal of Scientific Computing， 2018， 76（2）： 777.

［7］ 陳山. 求解波動方程的龍格—庫塔型方法及其地震波傳播模擬［D］. 北京. 清華大學， 2010.

［8］ HE Dongdong.An Unconditionally Stable Spatial Sixth-Order CCD-ADI Method for The Two-Dimensional Linear Telegraph Equation［J］. Numerical Algorithms， 2016， 72（4）： 1103.

［9］ MA Yuezhen， LI Xiangang， GE Y B.A High-accuracy Alternating Direction Implicit Method for Solving the Two-dimensional Wave Equation［J］. Journal of Sichuan Normal University： Natural Science， 2010， 33（2）： 179.

［10］魏劍英. 求解二維熱傳導方程的高精度緊致差分方法［J］. 西南師范大學學報：自然科學版， 2013， 38（12）： 5.

WEI Jianying.On a High-Order Compact Difference Method for Solving the Two-Dimensional Heat Equation［J］. Journal of Southwest University： Natural Science， 2013，38（12）： 5.

［11］YANG Dinghui， WANG Nan， CHEN Shan， at all. An Explicit Method Based on The Implicit Runge-Kutta Algorithm for Solving Wave Equations［J］. Bulletin of the Seismological Society of America， 2009， 99（6）： 3340.

［12］DAS S， LIAO Wenyuan， GUPTA A. An Efficient Fourth-Order Low Dispersive Finite Difference Scheme for A 2-D Acoustic Wave Equation［J］. Journal of Computational and Applied Mathematics， 2014， 258： 151.

［13］LI Jichun， CHEN Yitung， LIU Guoqing. High-Order Compact ADI Methods for Parabolic Equations［J］. Computers & Mathematics with Applications， 2006， 52（8/9）： 1343.

［14］MONHANTY R K. An Unconditionally Stable Finite Difference Formula for A Linear Second Order One Space Dimensional Hyperbolic Equation with Variable Coefficients［J］. Applied Mathematics and Computation， 2004， 165（1）： 229.

［15］WU Fengyan， CHENG Xiujun， LI Dongfang， et al. A Two-Level Linearized Compact ADI Scheme for Two-Dimensional Nonlinear Reaction-Diffusion Equations［J］. Computers & Mathematics with Applications， 2018： 2835.

［16］LIAO Wenyuan. On The Dispersion， Stability and Accuracy of A Compact Higher-Order Finite Difference Scheme for 3D Acoustic Wave Equation［J］. Journal of Computational and Applied Mathematics， 2014， 270： 571.

［17］ZHANG Wensheng， LI Tong， CHUNG E T. A New High Accuracy Locally One-Dimensional Scheme for The Wave Equation［J］. Journal of Computational and Applied Mathematics， 2011， 236（6）： 1343.

［18］GE Yongbin， CAI Zhiquan. High-order Compact Difference Schemes and Adaptive Method for Singular Degenerate Diffusion-Reaction Equations［J］. Chinese Journal of Computational Physics， 2017， 34（3）： 309.

［19］SAADANTMANDI A， DEHGHAN M. NumericalSolution of The One-Dimensional Wave Equation with An Integral Condition［J］. Numerical Methods for Partial Differential Equations， 2010， 23（2）： 282.

［20］BIAZAR J， GHAZVINI H.Homotopy Perturbation Method for Solving Hyperbolic Partial Differential Equations［J］. Computers & Mathematics with Applications， 2008， 56（2）： 453.

［21］姜蘊芝. 求解波動方程的高精度緊致顯式差分格式［D］.銀川：寧夏大學， 2016.

（編輯：溫澤宇）