# A High Accuracy Local One-Dimensional Explicit Compact Scheme for the 2D Acoustic Wave Equation

Authors: Mengling Wu; Yunzhi Jiang; Yongbin Ge
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

In this paper, we develop a highly accurate and efficient finite difference scheme for solving the two-dimensional (2D) wave equation. Based on the local one-dimensional (LOD) method and Padé difference approximation, a fourth-order accuracy explicit compact difference scheme is proposed. Then, the Fourier analysis method is used to analyze the stability of the scheme, which shows that the new scheme is conditionally stable and the Courant-Friedrichs-Lewy (CFL) condition is superior to most existing methods of equivalent order of accuracy in the literature. Finally, numerical experiments demonstrate the high accuracy, stability, and efficiency of the proposed method.


---

## Body

## 1. Introduction

Wave equation is mainly used to describe various wave phenomena in nature, such as sound wave, light wave, and water wave. Due to complex physical background, it is very difficult to obtain the exact solutions of practical problems. Therefore, researches on numerical solutions of initial value (or initial-boundary value) problems for wave equation have extremely important theoretical value and practical significance.The finite difference method is one of the most important and popular methods for solving the wave equation. It has many advantages, such as it is easy to be implemented, has low memory requirement, and has a high computing speed [1–8]. The second-order central difference scheme is the most commonly used, but when the number of spatial sampling points per wavelength is too small, the method has serious dispersion. In order to effectively suppress numerical dispersion, Yang et al. [9] proposed a local interpolation strategy based on the Runge-Kutta method. Later, for solving large-scale problems, Yang et al. [10] proposed an approximate analytic central difference method. However, both methods are conditionally stable; their maximum values of CFL conditions are 0.5080 [9] and 0.8440 [10], respectively. Most of the existing schemes with fourth-order accuracy in both space and time have very strict stability conditions. For instance, Wang et al. [11] developed a new finite difference stencil where the range of the CFL number is 0.7071. Feo et al. [12] combined novel tensor mimetic discretizations in space and a leapfrog approximation in time to produce an explicit scheme; the range of the CFL number is 0.8660. In addition, many implicit finite difference schemes are also used to solve the 2D wave equations [13–16]. For instance, Wang et al. [13] introduced two efficient fourth-order implicit time-space-domain finite difference schemes by using linear and nonlinear optimization methods. The high-order compact (HOC), compact Padé difference, and noncompact Padé difference implicit schemes were proposed in Ref. [14], all of which have a truncation error of Oτ4+τ2h2+h4. Since these implicit difference methods require iterations at each time step, this greatly increases the computational cost. In order to improve the computational efficiency, a combination of the finite difference method and splitting method is used by many researchers [17, 18]. Two types of classical splitting methods are the alternating direction implicit (ADI) method and the LOD method. Both of these methods simplify the multidimensional problems into a series of one-dimensional (1D) problems that form the tridiagonal systems and are easy to be solved. The ADI method was originally introduced by Peaceman and Rachford [19] to solve parabolic and elliptic equations. It was also used to solve hyperbolic equations [20, 21]. Up to now, many ADI methods have been studied to solve the wave equations [22–26]. For instance, Das et al. [22] and Liao et al. [23] used a combination of ADI and the Padé approximation to derive a difference scheme with fourth-order accuracy and low numerical dispersion. Their maximum values of CFL conditions are 0.7657 [22] and 0.7321 [23], respectively. Qin [24] constructed a compact Douglas scheme based on the ADI method, which has fourth-order accuracy in space, but only second-order accuracy in time.On the other hand, the LOD method was originally introduced by Samarskii [27] for solving the 2D hyperbolic equations. The LOD technique has been proven to be very useful in reducing computer memory and computational cost and improving stability [28, 29]. We notice that many researchers mainly study the homogeneous wave equation. For instance, Yun et al. [30] applied a new fourth-order accuracy optimal nearly analytic splitting method to solve the 2D acoustic wave equation based on the LOD method. Sim et al. [31] proposed a nearly analytic symplectic partitioned Runge-Kutta method based on the LOD technique, which has second-order temporal accuracy and fourth-order spatial accuracy. Very recently, Zhang [29] considered the nonhomogeneous wave equation and set an unknown parameter θ in front of the source term. Then, the value of θ was calculated by solving a series of large matrices, so as to obtain an implicit scheme with fourth-order accuracy in both time and space. This method increases computational complexity and reduces computational efficiency. Therefore, in this paper, we are aiming at developing an explicit compact difference method for solving the 2D nonhomogeneous wave equation based on the LOD method. The merits of the present method include high-order accuracy, better stability, and easy treatment for the nonhomogeneous term.The structure of this paper is arranged as follows. In Section2, the derivation of the new explicit HOC difference scheme based on the LOD method is given and the corresponding numerical algorithm is introduced. In Section 3, the stability of the new scheme is analyzed. In Section 4, the numerical examples are employed to verify our theoretical analysis results. Finally, the conclusion is given in Section 5.


## 2. Derivation of the New Scheme

The 2D wave equation with variable wave velocity is as follows:(1)∂2u∂t2=v2x,y∂2u∂x2+∂2u∂y2+fx,y,t,x,y,t∈Ω×0,T,with the initial conditions
(2)ux,y,0=φx,y,∂ux,y,0∂t=ψx,y,x,y∈Ω,and the Dirichlet boundary conditions
(3)ux,y,t=gx,y,t,x,y,t∈∂Ω×0,T,whereux,y,t represents pressure and vx,y is the wave velocity. fx,y,t is the source term. φx,y, ψx,y, and gx,y,t are known smooth functions. Ω=α,β×α,β, Ω⊂R2, and ∂Ω represent the boundary of Ω. α,β are constants, and α<β. 0,T is the time region.α,β is divided into N subintervals by equidistant grid, and the space step length is represented by h=β−α/N. 0,T is also equally divided into M subintervals, and the time step length is τ=T/M. We use xi,yj,tn which represents the mesh point, xi=α+ih,yj=α+jh,i,j=0,1,⋯N,tn=nτ,n=0,1,⋯M.The LOD method is used to ideally split the 2D wave equation (1) into two 1D wave equations (27) as follows:
(4)12∂2u∂t2=v2x,y∂2u∂x2+12fx,y,t,(5)12∂2u∂t2=v2x,y∂2u∂y2+12fx,y,t.To advance the solution fromuxi,yj,tn to uxi,yj,tn+1, we assume that Equation (4) holds from uxi,yj,tn to uxi,yj,tn+1/2 and Equation (5) holds from uxi,yj,tn+1/2 to uxi,yj,tn+1.Firstly, we consider Equation (4) at the grid point xi,yj,tn, i.e.,
(6)12∂2u∂t2i,jn=vi,j2∂2u∂x2i,jn+12fi,jn.For the time second derivative∂2u/∂t2 on the left side of Equation (6), we adopt the following expression:
(7)∂2u∂t2i,jn=4ui,jn+1/2−2ui,jn+ui,jn−1/2τ2−τ248∂4u∂t4i,jn+Oτ4.According to the Equation (4), we have
(8)∂4u∂t4=∂2∂t2∂2u∂t2=2v2∂2∂x2∂2u∂t2+∂2f∂t2.We define the central difference operator for the spatial second derivative as follows:(9)δx2ui,jn=ui+1,jn−2ui,jn+ui−1,jnh2.Substituting Equations (8) and (9) into Equation (7), we obtain
(10)∂2u∂t2i,jn=4ui,jn+1/2−2ui,jn+ui,jn−1/2τ2−vi,j2τ224∂2∂x2∂2u∂t2i,jn−τ248∂2f∂t2i,jn+Oτ4=4ui,jn+1/2−2ui,jn+ui,jn−1/2τ2−vi,j2τ224δx2∂2u∂t2i,jn−τ248∂2f∂t2i,jn+Oτ4+τ2h2=4ui,jn+1/2−2ui,jn+ui,jn−1/2τ2−vi,j2τ224h2∂2u∂t2i+1,jn+∂2u∂t2i−1,jn+vi,j2τ212h2∂2u∂t2i,jn−τ248∂2f∂t2i,jn+Oτ4+τ2h2.It can be rewritten as(11)ui,jn+1/2=2ui,jn−ui,jn−1/2+vi,j2τ496h2∂2u∂t2i+1,jn−2∂2u∂t2i,jn+∂2u∂t2i−1,jn+τ24∂2u∂t2i,jn+τ4192∂2f∂t2i,jn+Oτ4+τ2h2.Defining the grid ratioλ=τ/h, substituting Equation (4) into Equation (11), and omitting the truncation error term, we obtain
(12)ui,jn+1/2=2ui,jn−ui,jn−1/2+vi,j2τ2λ248vi+1,j2∂2u∂x2i+1,jn+vi−1,j2∂2u∂x2i−1,jn+vi,j2τ22−vi,j4τ2λ224∂2u∂x2i,jn+vi,j2τ2λ296fi+1,jn+fi−1,jn+τ24−vi,j2τ2λ248fi,jn+τ4192∂2f∂t2i,jn.The functionf and its derivative ∂2f/∂t2 are known at each grid point. Next, we use the fourth-order Padé scheme given in Ref. [32] to calculate the unknown term ∂2u/∂x2n in Equation (12) as follows:
(13)∂2u∂x2i+1,jn+10∂2u∂x2i,jn+∂2u∂x2i−1,jn=12ui+1,jn−2ui,jn+ui−1,jnh2+Oh4,i=1,2,⋯,N−1,j=0,1,⋯N,in which the nodes of the∂2u/∂x2i,jn on the boundaries can be obtained by using Equations (3) and (6):
(14)∂2u∂x20,jn=12v0,j2∂2u∂t2−f0,jn=12v0,j2∂2g∂t2−f0,jn,j=0,1,⋯N,(15)∂2u∂x2N,jn=12vN,j2∂2u∂t2−fN,jn=12vN,j2∂2g∂t2−fN,jn,j=0,1,⋯N.Then, we consider Equation (5) at the grid point xi,yj,tn+1/2, i.e.,
(16)12∂2u∂t2i,jn+1/2=vi,j2∂2u∂y2i,jn+1/2+12fi,jn+1/2.Using a similar treatment process as for Equation (6), omitting the truncation error terms, we can obtain
(17)ui,jn+1=2ui,jn+1/2−ui,jn+vi,j2τ2λ248vi,j+12∂2u∂y2i,j+1n+1/2+vi,j−12∂2u∂y2i,j−1n+1/2+vi,j2τ22−vi,j4τ2λ224∂2u∂y2i,jn+1/2+vi,j2τ2λ296fi,j+1n+1/2+fi,j−1n+1/2+τ24−vi,j2τ2λ248fi,jn+1/2+τ4192∂2f∂t2i,jn+1/2.The internal nodes of the∂2u/∂y2i,jn+1/2 can also be computed by the fourth-order Padé scheme [32].
(18)∂2u∂y2i,j+1n+1/2+10∂2u∂y2i,jn+1/2+∂2u∂y2i,j−1n+1/2=12ui,j+1n+1/2−2ui,jn+1/2+ui,j−1n+1/2h2+Oh4,i=0,1,⋯N,j=1,2,⋯,N−1,in which, the nodes of the∂2u/∂y2i,jn+1/2 on the boundaries can be obtained by using Equations (3) and (16):
(19)∂2u∂y2i,0n+1/2=12vi,02∂2u∂t2−fi,0n+1/2=12vi,02∂2g∂t2−fi,0n+1/2,i=0,1,⋯N,(20)∂2u∂y2i,Nn+1/2=12vi,N2∂2u∂t2−fi,Nn+1/2=12vi,N2∂2g∂t2−fi,Nn+1/2,i=0,1,⋯N.Applying Equations (12)–(20), we can complete the entire calculation from uxi,yj,tn to uxi,yj,tn+1 for Equation (1). The scheme is an explicit HOC difference scheme based on the LOD method, which can be denoted as the EHOC-LOD scheme. According to the derivation process, we know that the EHOC-LOD scheme is a four-level scheme with fourth-order accuracy in both time and space. Therefore, two start-up time steps need to be calculated at t=τ/2 and t=τ. To approximate the start-up time steps ui,j1/2 and ui,j1, using Taylor’s expansions, we have
(21)ui,j1/2=ui,j0+τ2∂u∂ti,j0+τ28∂2u∂t2i,j0+τ348∂3u∂t3i,j0+τ4384∂4u∂t4i,j0+Oτ5,(22)ui,j1=ui,j0+τ∂u∂ti,j0+τ22∂2u∂t2i,j0+τ36∂3u∂t3i,j0+τ424∂4u∂t4i,j0+Oτ5.For simplicity, we denote the time-independent functions as the following notations:(23)∂ς∂x=ςx,∂ς∂y=ςy,∂2ς∂x2=ςxx,∂2ς∂y2=ςyy,∂3ς∂x3=ςxxx,∂3ς∂x∂y2=ςxyy,∂3ς∂y3=ςyyy,∂3ς∂y∂x2=ςxxy,∂4ς∂x4=ςxxxx,∂4ς∂y4=ςyyyy,∂4ς∂x2∂y2=ςxxyy,in whichς can represent φ,ψ,v2, etc.By using Equations (1) and (2), we can get
(24)∂2u∂t2i,j0=vi,j2∂2u∂x2+∂2u∂y2i,j0+fi,j0=vi,j2φxx+φyyi,j+fi,j0,(25)∂3u∂t3i,j0=∂∂t∂2u∂t2i,j0=vi,j2∂2∂x2+∂2∂y2∂u∂ti,j0+∂f∂ti,j0=vi,j2ψxx+ψyyi,j+∂f∂ti,j0,(26)∂4u∂t4i,j0=∂2∂t2∂2u∂t2i,j0=∂2∂t2v2∂2u∂x2+∂2u∂y2+fi,j0=vi,j2∂2v2∂x2+∂2v2∂y2i,j∂2u∂x2+∂2u∂y2i,j0+2vi,j2∂v2∂xi,j∂3u∂x3+∂3u∂x∂y2i,j0+∂v2∂yi,j∂3u∂y3+∂3u∂y∂x2i,j0+vi,j42∂4u∂x2∂y2i,j0+∂4u∂x4+∂4u∂y4i,j0+vi,j2∂2f∂x2+∂2f∂y2i,j0+∂2f∂t2i,j0=vi,j2vxx2+vyy2i,jφxx+φyyi,j+2vx2i,jφxxx+φxyyi,j+2vy2i,jφyyy+φxxyi,j+vi,j4φxxxx+2φxxyy+φyyyyi,j+vi,j2∂2f∂x2+∂2f∂y2i,j0+∂2f∂t2i,j0.Substituting Equations (24)–(26) into Equations (21) and (22), respectively, and omitting the truncation error, we obtain
(27)ui,j1/2=φi,j+τ2ψi,j+τ28vi,j2φxx+φyyi,j+fi,j0+τ348vi,j2ψxx+ψyyi,j+∂f∂ti,j0+τ4384vi,j2vxx2+vyy2i,jφxx+φyyi,j+vi,j4φxxxx+φyyyyi,j+vi,j2∂2f∂x2+∂2f∂y2i,j0+∂2f∂t2i,j0+τ4192vi,j2vx2i,jφxxx+φxyyi,j+vy2i,jφyyy+φxxyi,j+vi,j2φxxyyi,j,(28)ui,j1=φi,j+τψi,j+τ22vi,j2φxx+φyyi,j+fi,j0+τ36vi,j2ψxx+ψyyi,j+∂f∂ti,j0+τ424vi,j2vxx2+vyy2i,jφxx+φyyi,j+vi,j4φxxxx+φyyyyi,j+vi,j2∂2f∂x2+∂2f∂y2i,j0+∂2f∂t2i,j0+τ412vi,j2vx2i,jφxxx+φxyyi,j+vy2i,jφyyy+φxxyi,j+vi,j2φxxyyi,j.Based on the description above, an algorithm for solving Equation (1) by using the EHOC-LOD scheme is given as follows:Step 1.
Letn=1, give the initial time step ui,j0, and the start-up time steps ui,j1/2 and ui,j1 are computed by using Equations (27) and (28) i,j=0,1,⋯N.Step 2.
The values of∂2u/∂x2i,jn are computed by Equations (13)–(15) i,j=0,1,⋯N. Then, the function values ui,jn+1/2 are computed by using Equation (12) i=1,⋯N−1,j=0,1,⋯N. The function values u0,jn+1/2 and uN,jn+1/2 are computed by Equation (3) j=0,1,⋯N.Step 3.
The values of∂2u/∂y2i,jn+1/2 are computed by Equations (18)–(20) i,j=0,1,⋯N. Then, the function values ui,jn+1 are computed by using Equation (17) i=0,⋯N,j=1,⋯N−1. The function values ui,0n+1 and ui,Nn+1 are computed by Equation (3) i=0,1,⋯N.Step 4.
Letn⟵n+1, repeat Steps 2 and 3 until time reaches the final moment, and the computation is terminated.


## 3. Stability Analysis

In this section, we analyze the stability of the EHOC-LOD scheme by using the Fourier analysis method.Lemma 1 (see [33]).
The sufficient and necessary condition for the roots of the quadratic equationμ2−bμ−c=0 with real coefficients to be not greater than one is c≤1, b≤1−c.Theorem 1.
The EHOC-LOD scheme in Section2 is stable if
(29)max1≤i,j≤Nvi,j·τh=vmaxλ≤0.8740,in which,vmax=max1≤i,j≤Nvi,j.Proof.
Letui,jn=ξneIσ1xieIσ2yj, uxxi,jn=ηneIσ1xieIσ2yj, and uyyi,jn=γneIσ1xieIσ2yj, where ξn,ηn,γn are amplitudes, σ1,σ2 are phase angles, and I=−1. We assume that fx,y,t is exact and does not produce error. Equations (13) and (18) can conclude that
(30)ηneIσ1xieIσ2yjeIσ1h+10+e−Iσ1h=12h2ξneIσ1xieIσ2yjeIσ1h−2+e−Iσ1h,(31)γn+1/2eIσ1xieIσ2yjeIσ2h+10+e−Iσ2h=12h2ξn+1/2eIσ1xieIσ2yjeIσ2h−2+e−Iσ2h.By using Euler’s formula,eIσh=cosσh+Isinσh,e−Iσh=cosσh−Isinσh, Equations (30) and (31) turn to
(32)ηnξn=12cosσ1h−1h2cosσ1h+5,(33)γn+1/2ξn+1/2=12cosσ2h−1h2cosσ2h+5.Letvi,jn+1/2=ui,jn and vi,j2max=a, and Equation (12) is written in matrix form:
(34)ui,jn+1/2vi,jn+1/2=2−110ui,jnvi,jn+aτ22−a2τ2λ224000uxxi,jnvxxi,jn+a2τ2λ248000uxxi+1,jnvxxi+1,jn+a2τ2λ248000uxxi−1,jnvxxi−1,jn.Substitutingui,jn,vi,jnT=ξneIσ1xieIσ2yj,uxxi,jn,vxxi,jnT=ηneIσ1xieIσ2yj into the formula above, we get
(35)ξn+1/2=2−110ξn+aτ22−a2τ2λ224000ηn+a2τ2λ248eIσ1h+e−Iσ1h000ηn.Substituting Equation (32) into Equation (35), we have
(36)ξn+1/2=2−110ξn+aτ22+a2τ2λ2cosσ1h−12400012cosσ1h−1h2cosσ1h+5ξn.Similarly, Equation (17) can be treated as
(37)ξn+1=2−110ξn+1/2+aτ22+a2τ2λ2cosσ2h−12400012cosσ2h−1h2cosσ2h+5ξn+1/2.Substituting Equation (36) into Equation (37), the error propagation matrix can be obtained as
(38)G=ξn+1ξn=B−110A−110=BA−1−BA−1,where
(39)A=2+12aλ2cosσ1h−1+a2λ4cosσ1h−122cosσ1h+5,B=2+12aλ2cosσ2h−1+a2λ4cosσ2h−122cosσ2h+5.The characteristic equation can be obtained as(40)μI−G=μ−BA−1B−Aμ+1=μ2−AB−2μ+1=0.According to Lemma1, b=AB−2≤2. Due to A and B having the same range of values, we only need to solve 0≤A≤2,0≤B≤2, or −2≤A≤0,−2≤B≤0.Case 1.
When0≤A≤2,0≤B≤2, we only need to solve 0≤A≤2, i.e.,
(41)−4cosσ1h+5≤12aλ2cosσ1h−1+a2λ4cosσ1h−12≤0.Firstly, for the right side of the inequality (41)
(42)12aλ2cosσ1h−1+a2λ4cosσ1h−12≤0,we have
(43)aλ2≤6.Then, for the left side of the inequality (41),
(44)−4cosσ1h+5≤12aλ2cosσ1h−1+a2λ4cosσ1h−12.Letting cosσ1h=ω,ω∈−1,1, the inequality (44) can be expressed as
(45)a2λ4ω2+12aλ2−2a2λ4+4ω+a2λ4−12aλ2+20≥0.LettingFω=a2λ4ω2+12aλ2−2a2λ4+4ω+a2λ4−12aλ2+20, we assume that Fω is a quadratic function of ω. Then, the function image is a parabola with an upward opening, the axis of symmetry is −12aλ2−2a2λ4+4/2a2λ4, the ordinate of the vertex is −12a2λ4−24aλ2−4/a2λ4. In this case, we only need to determine the position of the symmetry axis, so that Fmin≥0 on ω∈−1,1.Firstly, sinceaλ2≥0, the axis of symmetry −12aλ2−2a2λ4+4/2a2λ4≤1 always holds.Secondly, when the symmetry axis−12aλ2−2a2λ4+4/2a2λ4≤−1, we have 0≤aλ2≤3+13/2. Fω is an increasing function; we only need Fmin=F−1≥0, so aλ2≤3−5 or aλ2≥3+5. But 0≤aλ2≤3+13/2, so 0≤aλ2≤3−5≈0.7639.Finally, when the axis of symmetry−1≤−12aλ2−2a2λ4+4/2a2λ4≤1, we have aλ2≥3+13/2. We only need Fmin=F−12aλ2−2a2λ4+4/2a2λ4≥0, so −3−3/3≤aλ2≤−3+3/3. But aλ2≥3+13/2≈3.3028, so there is no solution.In summary, when0≤aλ2≤3−5, i.e., 0≤aλ≤3−5≈0.8740, the EHOC-LOD scheme is stable.Case 2.
When−2≤A≤0,−2≤B≤0, we only need to solve −2≤A≤0(46)−4≤12aλ2cosσ1h−1+a2λ4cosσ1h−122cosσ1h+5≤−2.Firstly, for the right side of the inequality (46),
(47)12aλ2cosσ1h−1+a2λ4cosσ1h−12≤−4cosσ1h+5.According the result of Case 1, we can get that
(48)aλ≥0.8740.Then, for the left side of the inequality (46),
(49)−8cosσ1h+5≤12aλ2cosσ1h−1+a2λ4cosσ1h−12.Letting cosσ1h=ω,ω∈−1,1, the inequality (49) be expressed as
(50)a2λ4ω2+12aλ2−2a2λ4+8ω+a2λ4−12aλ2+40≥0.LettingGω=a2λ4ω2+12aλ2−2a2λ4+8ω+a2λ4−12aλ2+40, we assume that Gω is a quadratic function of ω. Then, the function image is a parabola with an upward opening, the axis of symmetry is −12aλ2−2a2λ4+8/2a2λ4, and the ordinate of the vertex is 48a2λ4−192aλ2−64/4a2λ4. In this case, we only need to determine the position of the symmetry axis, so that Gmin≥0 on ω∈−1,1.Firstly, sinceaλ2≥0, the axis of symmetry −12aλ2−2a2λ4+8/2a2λ4≤1 always holds.Secondly, when the axis of symmetry−12aλ2−2a2λ4+8/2a2λ4≤−1, we have 0≤aλ2≤3+17/2. Gω is an increasing function; we only need Gmin=G−1≥0, so 0≤aλ2≤2 or aλ2≥4. But 0≤aλ2≤3+17/2≈3.5615, so 0≤aλ2≤2.Finally, when the axis of symmetry−1≤−12aλ2−2a2λ4+8/2a2λ4≤1, we have aλ2≥3+17/2. We only need Gmin=G−12aλ2−2a2λ4+8/2a2λ4≥0, so aλ2≥6+43/3. But aλ2≥3+17/2, so aλ2≥3+17/2≈3.5615.In summary, due toaλ≥0.8740∩aλ≤2∩aλ≥3+17/2=∅, there is no solution for this case.Combining Cases1 and 2, we can conclude that when aλ=vmaxλ≤0.8740, the EHOC-LOD scheme is stable. This completes the proof.Table1 gives the stability conditions of various schemes in Refs. [9–12, 15, 16, 22, 23, 28, 30] and the present EHOC-LOD scheme for the 2D wave equation. Among them, some methods are aimed at the constant wave number wave equation such as Refs. [9, 10, 12, 16, 22, 28, 30]; v represents the absolute value of the constant wave coefficient. The others deal with the wave equation with variable wave coefficients such as Refs. [11, 15, 23]; vmax=maxα≤x,y≤βvx,y represents the maximum value of the variable wave coefficients. The theoretical analysis result of the stability condition of the EHOC-LOD scheme in this paper is vmaxλ∈0,0.8740. So, it is more flexible than all the other schemes listed in the literature.Table 1
The stability conditions of various finite difference schemes for solving the 2D wave equations.
SchemesStability conditionNACD [10]vλ∈0,0.5080ONAS [30]vλ∈0,0.5629FD [16]vλ∈0,0.7071LOD [28]vλ∈0,0.7321NCPD-ADI [22]vλ∈0,0.8186IRK-DSM [9]vλ∈0,0.8440Mimetic [12]vλ∈0,0.8660FD [15]vmaxλ∈0,0.6124N-TSD-TEM [11]vmaxλ∈0,0.7070NCV-CPD-ADI [23]vmaxλ∈0,0.7321EHOC-LODvmaxλ∈0,0.8740


## 4. Numerical Experiments

This section uses some numerical examples to confirm the accuracy and stability of the proposed EHOC-LOD method. All programs are written in Fortran 90 language by using double-precision arithmetic and run on an Intel Core i5-4210U CPU@1.70 GHz 2.40 GHz desk computer with 4 GB memory.L2 and L∞ norm errors and the rate of convergence are obtained by the following definitions:
(51)L∞=max0≤i,j≤Nui,jM−uxi,yj,tM,L2=h2∑i,j=0Nui,jM−uxi,yj,tM2,Rate=logL∞h1/L∞h2logh1/h2.The CPU time for all the computation is given in seconds.Example 1 (see [28]).
(52)∂2u∂t2=∂2u∂x2+∂2u∂y2,0≤x,y≤1,t>0,ux,y,0=sinπxsinπy,∂ux,y,0∂t=0,u0,y,t=0,u1,y,t=0,ux,0,t=0,ux,1,t=0.
The exact solution is(53)ux,y,t=cos2πtsinπxsinπy.Example1 is a 2D homogeneous wave equation with constant coefficients. Table 2 shows the L2 and L∞ norm errors for the EHOC-LOD scheme with τ=0.0002 and T=0.02. In order to compare, we give the calculated results of the LOD scheme proposed by Zhang et al. in Ref. [28] in Table 2. The results show that the L∞ and L2 norm errors of the EHOC-LOD scheme are smaller than those of the LOD scheme [28]. In order to verify the computed accuracy of the EHOC-LOD scheme in space direction, in Figure 1, we draw log-log plots for the L2 and L∞ norm errors in the space direction for the LOD scheme [28] and the EHOC-LOD scheme. That displays that the slope of the line of the EHOC-LOD scheme is close to 4, which means that the EHOC-LOD scheme can converge with the fourth order in the space direction. When the CFL number is fixed to 0.2 and T=0.02, the L2 and L∞ norm errors are listed in Table 3. Obviously, the numerical solution of the EHOC-LOD scheme is a bit better than that of the LOD scheme [28]. Then, two log-log plots for the L2 and L∞ norm errors in the time direction for the LOD scheme [28] and the EHOC-LOD scheme are shown in Figure 2. Figure 2 shows that the slope of the line of the EHOC-LOD scheme exceeds 4, which means that the EHOC-LOD scheme can converge with the fourth order in the time direction.Example 2 (see [24]).
(54)∂2u∂t2=∂2u∂x2+∂2u∂y2+3π2e−πtsinπx+πy,ux,y,0=sinπx+πy,∂ux,y,0∂t=−πsinπx+πy,u0,y,t=e−πtsinπy,u1,y,t=e−πtsinπ+πy,ux,0,t=e−πtsinπx,ux,1,t=e−πtsinπx+π.
The exact solution is(55)ux,y,t=e−πtsinπx+πy.Table 2
L2 and L∞ norm errors when T=0.02 and τ=0.0002 for Example 1.
hLOD scheme [28]EHOC-LOD schemeL∞L2L∞L21/103.8515(-6)1.7507(-6)1.5824(-7)7.9120(-8)1/202.4002(-7)1.1429(-7)9.8614(-9)4.9307(-9)1/401.4990(-8)7.3122(-9)6.1587(-10)3.0792(-10)1/809.3644(-10)4.6244(-10)3.8498(-11)1.9238(-11)1/1605.8505(-11)2.9068(-11)2.4473(-12)1.2005(-12)1/3204.7895(-12)2.3824(-12)1.8097(-13)7.5502(-14)Figure 1
Log-log plots for theL2 (a) and L∞ (b) norm errors in space direction.
(a)(b)Table 3
L2 and L∞ norm errors when T=0.02 and the CFL number is fixed at 0.2 for Example 1.
τLOD scheme [28]EHOC-LOD schemeL∞L2L∞L20.00405.5715(-8)2.7858(-8)5.3806(-8)2.6903(-8)0.00254.8646(-9)2.4323(-9)4.6971(-9)2.3486(-9)0.00201.3816(-9)6.9080(-10)1.3333(-9)6.6664(-10)0.00163.8135(-10)1.9071(-10)3.6722(-10)1.8363(-10)0.00102.5111(-11)1.2545(-11)2.3096(-11)1.1537(-11)0.00087.5682(-12)3.7770(-12)6.1298(-12)3.0595(-12)Figure 2
Log-log plots for theL2 (a) and L∞ (b) norm errors in time direction.
(a)(b)Example2 is a 2D nonhomogeneous wave equation with constant coefficients. In order to make a better comparison, we program the algorithms of the CPR, CDOU, and CSchemeI schemes with second-order accuracy in time and fourth-order accuracy in space in Ref. [24] also in Fortran 90 language. We run the programs on the same computer to make the comparison fair. When T=2, for different spatial grid step length, in Table 4, the L2 norm errors and CPU time are shown by using the present EHOC-LOD scheme with τ=0.5h, the CPR, CDOU, and CSchemeI schemes with τ=h2 in Ref. [24]. We can notice that the EHOC-LOD scheme obtains fourth-order accuracy in both time and space, but those schemes in Ref. [24] must take τ=Oh2 to achieve fourth-order accuracy in the space direction. Moreover, the present method gets a more accurate solution. When we define the time as T and the number of grids as N, the total number of time advance steps of the method in this paper is 2NT, while the total number of time advance steps of those methods in Ref. [24] is N2T; therefore, the EHOC-LOD scheme expends much less CPU time than those methods in Ref. [24]. With the increase of time T and grid number N, the method in this paper has more significant advantages.Table 4
TheL2 norm error and CPU time at T=2 for Example 2.
SchemesτN=40N=60N=120L2CPUL2CPUL2CPUCPR [24]h23.3075(-7)2.39906.5916(-8)11.89201.5789(-9)200.8900CDOU [24]h23.3101(-7)3.56606.5886(-8)15.69601.9432(-9)218.1560CSchemeI [24]h21.0870(-6)3.60502.1494(-7)17.55108.2722(-8)270.0030EHOC-LOD0.5h2.0132(-8)0.52404.1555(-9)1.65602.7109(-10)13.3130Then, we test the stability condition of the EHOC-LOD method. TheL∞ norm error of the EHOC-LOD scheme with h=1/32 for different T and τ is shown in Table 5. As the time T increases, when τ≤1/37,vmaxλ≤0.8649, the EHOC-LOD scheme is convergent. When τ≥1/36,vmaxλ≥0.8889, the EHOC-LOD scheme is divergent. These computed results are in good agreement with the theoretical analysis.Example 3 (see [23]).
(56)∂2u∂t2=1+sin2x+sin2y∂2u∂x2+∂2u∂y2+3+2sin2x+2sin2ye−tcosxcosy,0≤x,y≤π,t>0,ux,y,0=cosxcosy,∂ux,y,0∂t=−cosxcosy,u0,y,t=e−tcosy,uπ,y,t=−e−tcosy,ux,0,t=e−tcosx,ux,π,t=−e−tcosx.
The exact solution is(57)ux,y,t=e−tcosxcosy.Table 5
TheL∞ norm error with h=1/32 for different τ and T for Example 2.
τvmaxλT=1T=5T=101/550.58182.3265(-8)5.9423(-8)4.9813(-8)1/500.64001.4417(-8)5.6307(-8)1.0922(-7)1/450.71114.6947(-8)7.4385(-8)1.6359(-7)1/400.80001.1738(-7)1.1719(-7)3.2546(-7)1/380.84211.5834(-7)1.6067(-7)3.8232(-7)1/370.86491.8287(-7)1.6693(-7)4.1226(-7)1/360.88891.7126(-5)3.8468+0118.8343+0331/350.91433.7430(-2)5.1554+0294.3232+0681/340.94123.71291.8085+0407.3187+0891/330.9697154.65905.2153+0487.3864+106WhenT=1, for various h and τ, we compute the L∞ norm error and Rate by using the NCV-CPD-ADI scheme in Ref. [23] and the EHOC-LOD scheme in Table 6. The EHOC-LOD scheme produces the numerical solution with a smaller error than the NCV-CPD-ADI scheme does [23]. In addition, The EHOC-LOD scheme can achieve the fourth-order accuracy in both time and space, while the accuracy of the NCV-CPD-ADI scheme [23] is lower than that of the fourth order.Table 6
L∞ norm error and Rate when T=1 with various τ and h for Example 3.
h,τNCV-CPD-ADI [23]EHOC-LODL∞RateL∞Rateπ/40,1/802.7242(-6)2.5481(-8)π/80,1/1602.2523(-7)3.59641.5691(-9)4.0298π/160,1/3201.7543(-8)3.68249.7282(-11)4.0161π/320,1/6401.2928(-9)3.76236.0591(-12)4.0069The stability condition is given byvmaxλ≤0.8740, where vmax=maxx,y∈0,π×0,π1+sin2x+sin2y=3 for this example. Therefore, the grid ratio allowed by the EHOC-LOD scheme is τ/h≤0.8740/3≈0.5046. The grid ratio allowed by the NCV-CPD-ADI scheme [23] is τ/h<0.7321/3≈0.4226. When we fix h=π/200, then τ<0.4226π/200≈0.00663 and 1/150<0.00663<1/151. We take several numbers near 1/150; Table 7 shows the results. When τ≤1/149, the L∞ norm error of the NCV-CPD-ADI scheme [23] and the present scheme are both quite accurate. When τ≥1/148, the L∞ norm error of the NCV-CPD-ADI scheme [23] increases faster, and eventually, it is divergent due to exceeding the stability range. However, the L∞ norm error of the EHOC-LOD scheme is still quite accurate and convergent. Therefore, the CFL condition of the scheme in this paper is superior to that of the NCV-CPD-ADI scheme [23].Table 7
TheL∞ norm error with h=π/200 at T=1 for different τ for Example 3.
τNCV-CPD-ADI [23]EHOC-LOD1/1518.1293(-9)3.7842(-11)1/1508.1398(-9)3.7569(-11)1/1498.1505(-9)3.7639(-11)1/1481.2298(-7)3.7398(-11)1/1470.033473.7501(-11)To further verify the stability range of the EHOC-LOD scheme, Table8 shows the L∞ norm error of the EHOC-LOD scheme with h=π/200 for different T and τ. We notice that as time T grows, when τ≤1/127,vmaxλ≤0.8682, the EHOC-LOD scheme is convergent. When τ≥1/125,vmaxλ≥0.8821, the EHOC-LOD scheme is divergent. The numerical experiment results are consistent with the theoretical ones.Example 4 (see [23]).
Example 4 is a0,5km×0,5km two-layer model with a wave source. When y∈0,2500m, the vx,y is 1500m/s; when y∈2500m,5000m, the vx,y is 3000m/s. Ricker’s wavelet source that generates the wave is given by
(58)fx,y,t=δx−x0,y−y01−2π2fp2t−dr2e−π2fp2t−dr2,where δx−x0,y−y0 is the Dirac delta distribution, x0=2500m, and y0=2000m. fp=10Hz is the peak frequency. dr=0.5/fp is the temporal delay that is used to ensure zero initial conditions. The uniform grid size h=hx=hy=12.5m and Δt=0.001s are chosen to solve Example 4.Table 8
TheL∞ norm error with h=π/200 for different T and τ for Example 3.
τvmaxλT=1T=5T=101/1700.61263.9454(-11)2.9495(-11)7.5373(-12)1/1600.68923.8399(-11)2.8656(-11)6.6731(-12)1/1500.73513.7968(-11)2.8299(-11)6.4342(-12)1/1400.78763.6524(-11)2.7555(-11)6.7108(-12)1/1310.84173.4756(-11)2.6059(-11)7.1994(-12)1/1290.85483.5237(-11)2.5741(-11)7.4449(-12)1/1270.86823.3824(-11)2.5329(-11)7.5574(-12)1/1250.88213.6304(-11)8.3126+0352.9323+0911/1230.89652.54293.5497+075OverflowWe draw snapshots of the wave field at different times. Figure3 describes wave field snapshots computed by the EHOC-LOD scheme: (a) t=0.3s, (b) t=0.65s, (c) t=0.95s, and (d) t=1.15s. The distance from the source to the rock area is 2500m‐2000m=500m, so when t=500m/1500m/s=0.3333s, the wave will reach the rock area. In Figure 3(a), the wave still constitutes a complete circle until it reaches the interface of two layers. As time grows, the wave has reached the rock area; the reflection phenomenon is occurring near the interface, as shown in Figure 3(b). At this time, we can see that the upper and lower parts are two circles of different sizes, because the wave velocity of the two layers is constant. When the time is 0.95 s, the reflection phenomenon and refraction phenomenon will be further generated, forming a circular wave surface, as shown in Figure 3(c). Finally, in Figure 3(d), the left and right sides of the wave are slightly impacted, resulting in reflection.Figure 3
Wavefield snapshots: (a)t=0.3s; (b) t=0.65s; (c) t=0.95s; (d) t=1.15s.
(a)(b)(c)(d)


## 5. Conclusions

In this paper, a high accuracy explicit compact difference method for solving the 2D wave equation is presented. Firstly, we use the LOD method to transfer the 2D wave equation into two 1D problems. Then, combined with the fourth-order Padé method, an EHOC-LOD difference scheme is established for the two 1D problems with fourth-order accuracy in both time and space. Afterwards, the stability of the EHOC-LOD scheme is analyzed by the discrete Fourier method, and the range of the stability conditionvmaxλ∈0,0.8740 is obtained. Finally, numerical experiments are conducted to verify the accuracy and stability of the present method. The advantages of this method are mainly reflected in the following aspects:
(i)
It can be used not only to compute constant coefficient problems but also to compute variable coefficient problems, even if the wave velocity is discontinuous(ii)
It can be used not only to compute homogeneous problems but also to compute nonhomogeneous problems(iii)
Because of the use of the LOD technique, the computational complexity of the proposed method is significantly reduced. Furthermore, each 1D problem is solved by an explicit computational method, and the computational efficiency is further improved(iv)
The proposed method achieves the fourth-order accuracy; the stability range of the method is the best among all those schemes of equivalent order of accuracy in the literatureRecently, some authors consider HOC difference schemes for the 3D wave equations [34–36] and the Rosenau-Regularized Long-Wave (RLW) equation [37, 38]. The ADI method was used in Ref. [34, 36] and Mohanty and Gopal [35] proposed an implicit difference scheme. For the 3D Rosenau-RLW equation, Li [37] developed a fourth-order compact and energy-conservative scheme, which is a two-level nonlinear scheme. For the generalized Rosenau-RLW equation, Dimitrienko et al. [38] proposed a conservative fourth-order compact scheme, which is a three-level linear scheme. Generalizing the present EHOC-LOD method to these problems is our ongoing research work. The research results will be reported in the near future.



---

*Source: 9743699-2022-05-26.xml*