An explicit fourth-order compact difference scheme for solving the 2D wave equation

In this paper, an explicit fourth-order compact (EFOC) difference scheme is proposed for solving the two-dimensional(2D) wave equation. The truncation error of the EFOC scheme is O(τ 4 + τ 2h2 + h4), i.e., the scheme has an overall fourth-order accuracy in both time and space. Because the scheme is explicit, it does not need any iterative processes. Afterwards, the stability condition of the scheme is obtained by using the Fourier analysis method, which has a wider stability range than other explicit or alternation direction implicit (ADI) schemes. Finally, some numerical experiments are carried out to verify the accuracy and stability of the present scheme.


Introduction
In this paper, we consider the 2D wave equation as follows: with the initial conditions u(x, y, 0) = ϕ(x, y), (x, y) ∈ Ω, ( 2 ) ∂u(x, y, 0) and the boundary conditions u(0, y, t) = α 0 (y, t), u(l, y, t) = α l (y, t), (4) u(x, 0, t) = β 0 (x, t), u(x, l, t) = β l (x, t), (x, y, t) ∈ ∂Ω × [0, T], where Ω = (x, y) : 0 ≤ x, y ≤ l and ∂Ω is the boundary of Ω. Unknown function u(x, y, t) is displacement, a is the wave velocity, f (x, y, t) is the source term, ϕ(x, y) and ψ(x, y) are initial displacement and speed, respectively, α 0 (y, t), α l (y, t), β 0 (x, t), and β l (x, t) are displacements on the boundaries. We assume that all the functions are sufficiently smooth for achieving the accuracy order of the difference scheme. Solution to wave propagation problems has received considerable attention for many years [1-17, 19-32, 34]. Numerical methods are particularly attractive for structurally complex subsurface geometries because of the great difficulties encountered in obtaining analytical solutions. Among the numerical techniques available in the literature, the method of finite differences is particularly versatile [1-15, 19-23, 25-28, 31, 32, 34]. The most commonly used for solving this equation is the second-order central difference scheme. Since it is lower accuracy and lower resolution, at least 20 grid points is necessary to be input into each propagating wavelength. In other words, if fewer grid points are used, dispersion errors would make numerical solution oscillating. In order to cut down dispersion error and get dependable results, people usually have to use fine grids for computation, which inevitably increases computational costs and storage requirements significantly. An alternate approach for overcoming this difficulty is pseudospectral (PS) method [17]. In theory of the PS method, only two grid points are required for each wavelength (although four or more are required in practice), and discrete Fourier transforms on each time step are quite expensive. The space operators of the PS method fit the Nyquist frequency, but it needs Fourier transform that is time-consuming and the use of a global basis, resulting in inaccuracies in wave fields for models with strong heterogeneity or sharp boundaries [16]. To effectively eliminate numerical dispersion error when coarse grids are used to discretize the wave equation, Yang et al. [30] proposed the nearly-analytic discrete (NAD) approach. Later, they [29,31] further developed several NAD-type approaches to restrain numerical dispersion. Recently, many research works [15,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] have been done to develop higher-order finite difference methods for the above mentioned wave equations. Furthermore, it has been demonstrated that high accuracy numerical approaches and techniques are very effective in restraining the dispersion errors [5,8,21,31].
On the other hand, for multi-dimensional problems, operator splitting methods, such as the alternating direction implicit (ADI) [5,6,13,[19][20][21][22][23]25] methods and local onedimensional (LOD) approaches [15,34], are studied extensively. The advantage of the operator splitting methods is the solution of a multi-dimensional problem reduced into a series of 1D problems, and thus the efficient algorithms for 1D such as tri-diagonal matrix algorithm can be employed, which is more efficient than full implicit schemes. Unconditionally stable high-order compact ADI schemes with spatially fourth-order and temporally second-order are investigated in [6,13,22,23,25]. Liao and Sun [19,20] used Richardson extrapolation technique to promote the temporally second-order to fourthorder accuracy, but the second-order accuracy solutions on finer grid are needed to get the fourth-order solution on coarser grid. So, computed results on fine grid still show second-order accuracy. Das et al. [23] proposed some ADI schemes for solving the 2D wave equation, which are fourth-order accuracy in both time and space. However, they are conditionally stable and suffer from the stability conditions which only allow Courant-Friedrichs-Lewy (CFL) numbers from 0.7321 to 0.8186. Liao and Sun [19] extended this method to the 3D wave equation. It is also fourth-order accuracy in both time and space and the stability condition allows CFL number to be about 0.6079. For the study of LOD methods, the readers are referred to [15,34], which still suffers from a strict stability condition for the fourth-order scheme in both time and space.
In this paper, we shall introduce an explicit high-order compact difference scheme for solving the 2D wave equation. Our basic strategy is to perform a modified equation technique for temporal derivative, but we use the fourth mixed derivative of temporal and spatial variables instead of pure fourth derivative of spatial variables, thus only three grid points are used for spatial dimension. And spatial derivatives are computed by directly using the classical fourth-order Padé schemes [18] at known time steps. This allows us to obtain an explicit difference scheme that is formally fourth-order accurate in both space and time. An outline of the article is as follows. In the next section, a new three-level explicit high-order compact difference scheme for solving Eq. (1) is introduced. Section 3 gives the Fourier analysis about the scheme, and the range of stability condition is obtained. In Sect. 4, we give some numerical experiments to show numerical stability and accuracy of the present new scheme. The last section draws conclusions.

EFOC difference scheme
Mesh point is marked as (x i , y j , t n ), in which x i = ih, y j = jh, t n = nτ , i, j = 0, 1, 2, . . . , N , and n > 0. Let u n ij be an approximation to u ( x i , y j , t n ). τ is temporal mesh step length and h is spatial mesh step length.

Computation of the first time step
Using Taylor's expansions, we have Employing Eq. (1) and initial conditions (2) and (3), we have Substituting (7)-(9) into (6) and throwing away the truncation error term O(τ 5 ), we obtain the fourth-order compact difference scheme for the first time step as follows: For interior grid points at the (n)th time moment, we use the following fourth-order Padé approximation schemes [18]: and for grid points on the boundaries, we use Eqs. (1), (4), and (5) to give out directly We notice that the coefficient matrixes of Eqs. (11) and (12), with boundary conditions (13)- (16), are tri-diagonal, so they can be solved by the efficient tri-diagonal matrix algorithm.

Total time marching scheme
For the temporal derivative ( ∂ 2 u ∂t 2 ) n ij , we use the second-order central difference and keep the leading term of the truncation errors to get , and by using Eq. (1), we get Substituting (18) into (17), we have in which δ 2 . Rearranging (19) and neglecting the truncation error, we have Using Eq. (1) again, (20) becomes 4 12 Equation (21) is an explicit high-order compact difference scheme for solving Eq. (1). Through the derivation process, it is easy to know that the truncation error of the scheme is O(τ 4 + τ 2 h 2 + h 4 ), i.e., it has the fourth-order accuracy for both time and space. We mark it as EFOC scheme (explicit fourth-order compact scheme). It should be pointed out that ( ∂ 2 u ∂x 2 ) n ij and ( ∂ 2 u ∂y 2 ) n ij must be computed by Eqs. (11)-(16) in advance. Now, we can get an algorithm for solving wave equation (1) with initial and boundary conditions (2)-(5) as follows.

Stability analysis
In this part, von Neumann linear stability analysis method is used to get the stability condition of the present EFOC scheme (21). We suppose that u is a periodic function, v n+1 ij = u n ij , and the source function f ≡ 0. Letting u n ij = ξ n e I(σ 1 x i +σ 2 y j ) , (u xx ) n ij = η n e I(σ 1 x i +σ 2 y j ) , (u yy ) n ij = γ n e I(σ 1 x i +σ 2 y j ) , in which ξ , η, and γ are amplitudes, σ 1 and σ 2 are wave numbers in x and y direction, respectively. I = √ -1 is an imaginary unit.

Numerical experiments
To validate the accuracy and stability of the present EFOC scheme, two test problems are solved by the present EFOC scheme with different temporal step length, spatial step length, and final time. Numerical results for L ∞ and L 2 norm errors and convergence rate are given, in which L ∞ and L 2 norm errors and convergence rate are defined as follows: .  the exact solution is

Problem 1 ([5])
. Table 2 and Table 3 give L ∞ and L 2 norm errors for t = 1 with τ = 0.0025 and different spatial grid sizes for Problem 1 with the EFOC scheme, respectively. For comparison, we also list the computed results in Ref. [5], in which the THOC-ADI scheme, the HOC-LOD scheme, the CPD-ADI scheme, and the IPD-ADI schemes were given. We find that all of these schemes reach fourth-order accuracy. The EFOC scheme yields a bit more accurate solution than the THOC-ADI scheme, the HOC-LOD scheme, and the CPD-ADI scheme, but a bit less accurate than the IPD-ADI scheme. But we notice that the IPD-ADI scheme, with five grid points along one direction, is non-compact, so the computation process is more complicated and the computational cost is higher than that of the present EFOC scheme. Figure 1 shows the numerical solution (a), the exact solution (b), the absolute error (c), and the contour plots of the numerical solution and exact solution (d) by the EFOC scheme in this article when τ = 0.0025, t = 1, N = 40, respectively, for Problem 1. It can be seen from Fig. 1 that the numerical solution in this article agrees well with the exact solution. u(x, y, t) = sin(πx) sin(πy) cos( √ 2πt).
We computed this problem with the present EFOC scheme and compared the results with Ref. [34], in which a local one-dimensional scheme (New LOD scheme) and a typical  Table 4 The L ∞ and L 2 norm errors with τ = 0.0002 at t = 0.02 for different h for Problem 2 h New LOD scheme [34] Typical 4th-order scheme [34] E F O Cs c h e m e L 2 -error L ∞ -error L 2 -error L ∞ -error L 2 -error L ∞ -error 1/10 1.75(-6) 3.85(-6) 2.32 (-6) fourth-order scheme are developed. Table 4 gives the numerical results with τ = 0.0002, h = 1/10, 1/20, . . . , 1/320, and t = 0.02. It shows that the L ∞ and L 2 norm errors of the three schemes decrease on behalf of the spatial grid size and fourth-order accuracy is obtained. We notice that the present EFOC scheme produces more accurate solution than the new LOD scheme and the typical fourth-order scheme [34]. Figure 2 shows the numerical solution (a), the exact solution (b), the absolute error (c), and the contour plots of the numerical solution and exact solution (d) for the EFOC scheme of Problem 2 when τ = 0.0002, t = 0.02, N = 40, respectively. In Fig. 2(d) the solid line is exact solution and the dotted line is numerical solution. It is not difficult to find that the numerical solution obtained from the EFOC scheme in this article agrees well with the exact solution, and the dispersion error is very small.