1 Introduction
Finite difference (FD) method has been widely used in various science and engineering applications for several reasons such as easy implementation, high efficiency etc, when the analytical solution is not available. One example is the acoustic wave equation when a nonzero point source function is included in the equation. In particular, the highorder FD methods have attracted the interests of many researchers working on seismic modelling (see Chen2007 ; Cohen1996 ; Liu2009a ; Takeuchi2000 and references therein) due to the highorder accuracy and effectiveness in suppressing numerical dispersion. Moreover, the number of grid points per wavelength required by higher order FD methods is significantly less than that of the conventional FD methods.
Recently, a great deal of efforts have been devoted to develop highorder FD schemes for the acoustic equations, and many accurate and efficient methods have been developed. Levander Levander1988 addressed the costeffectiveness of solving real problems using highorder spatial derivatives to allow a more coarse spatial sample rate. In Liu2009a , the authors used a plane wave theory and the Taylor series expansion to develop a low dispersion timespace domain FD scheme with error O for 1D, 2D and 3D acoustic wave equations, where and represent the time step and spatial grid size, respectively. It was then shown that, along certain fixed directions the error can be improved to . In Cohen1996 , Cohen and Poly extended the works of Dablain Dablain1986 , Shubin and Bell Shubin1987 and Bayliss et al. Bayliss1986 and developed a fourthorder accurate explicit scheme with error of O to solve the heterogeneous acoustic wave equation. Moreover, it has been reported that highly accurate numerical methods are very effective in suppressing the annoying numerical dispersion Finkelstein2007 . Highorder FD method is of particular importance for largescale 3D acoustic wave equation, as it requires less grid pointsEtgen2007 .
These methods are accurate but are noncompact, which give rise to two issues: efficiency and difficulty in boundary condition treatment. For example, the conventional noncompact fourthorder FD scheme requires a fivepoint stencil in 1D to approximate , while in 2D, a 9point stencil is needed in approximating . In 3D problem, a 13point stencil is required to approximate the derivative
with fourthorder accuracy. To overcome these difficulties, a variety of compact higherorder FD schemes have been developed for hyperbolic, parabolic and elliptical partial differential equations. In
Chu1998 , the authors developed a family of fourthorder threepoint combined difference schemes to approximate the first and secondorder spatial derivatives. In Kim2007 , the authors introduced a family of threelevel implicit FD schemes which incorporate the locally onedimensional method. For more recent compact higherorder difference methods, the readers are referred to Shukla2005 .For threedimensional problems, an implicit scheme results in a block tridiagonal system, which is solved at each time step. Direct solution of such large block linear system is very inefficient, therefore, some operator splitting techniques are used to convert the threedimensional problem into a sequence of onedimensional problems. One widely used method is the ADI method, which was originally introduced by Peaceman and Rachford Peaceman1955 to solve parabolic and elliptic equations. Since then a lot of developments have been made over the years for hyperbolic equations Douglas1966 ; Lees1962 . Later, Fairweather and Mitchell Fairweather1965 developed a fourthorder compact ADI scheme for solving the wave equation. Some other related work can be found in Li1991 ; Ristow1997 . Combined with Padé approximation of the finite difference operator, some efficient and highorder compact finite difference methods have been developed to solve the acoustic wave equations in 2D and 3D with constant velocity Das2014 ; Liao2014 . However when the velocity is a spatially varying function, it is difficult to apply this technique because the algebraic manipulation is not applicable here. Nevertheless, some research work on accurate and lowdispersion numerical simulation of acoustic wave propagation in heterogeneous media have been reportedZhang2011 , in which a layered model consisting of multiple horizontal homogeneous layers was considered, with the method development and stability analysis are based on constant velocity model.
In Liao2018 a new fourthorder compact ADI FD scheme was proposed to solve the twodimensional acoustic wave equation with spatially variable velocity. Here we extended this method and the Padé approximation based highorder compact FD scheme in Das2014 to the 3D acoustic wave equation with nonconstant velocity. The new method is compact and efficient, with fourthorder accuracy in both time and space. To our best knowledge, it is the first time that the energy method has been used to conduct theoretical stability analysis of the Padéapproximation based compact scheme. It has been shown that the obtained stability condition (CFL condition) is highly consistent to the CFL condition obtained by other methods. The rest of the paper is organized as the follows. We first give a brief introduction of the acoustic wave equation and several existing standard secondorder central FD schemes, compact highorder FD schemes for spatial derivatives and some other related highorder method in Section 2, then derive the new compact fourthorder ADI FD scheme in Section 3. Stability analysis of the new method is presented in Section 4, which is followed by three numerical examples in Section 5. Finally, the conclusions and some future works are discussed in Section 6.
2 Acoustic wave equation and existing algorithms
Consider the 3D acoustic wave equation
(1)  
(2)  
(3)  
(4) 
where represents the wave velocity. Here is a finite computational domain and is the source function. We denote for the sake of simple notation.
First assume that is a 3D rectangular box : , which is discretized into an grid with spatial grid sizes , and . Let be the time stepsize and denote the numerical solution at the grid point and time level . Let’s first define the standard central difference operators
The standard secondorder central FD schemes are then given by approximating the second derivatives in Eq.(1). For example,
(5)  
(6) 
The approximations of and are similar to Eq. (6).
To improve the method to fourthorder in space, the conventional highorder FD method was derived by approximating the derivatives using more than three points in one direction, which results in larger stencil. For instance, if points are used to approximate , one can obtain the following formula
(7) 
which can be as accurate as order in . The conventional highorder FD method is accurate in space but suffers severe numerical dispersion. Another issue is that it requires more computer memory due to the larger stencil for implicit method. Moreover, more points are needed to approximate the boundary condition.
To improve the accuracy in time, a class of timedomain highorder FD methods have been derived by Liu and Sen Liu2009a . The idea of the timedomain highorder FD method is to determine coefficients using timespace domain dispersion. As a result, the coefficient will be a function of . It was noted that in 1D case, the timedomain highorder FD method can be as accurate as order in both time and space, provided some conditions are satisfied, while for multidimensional case, order is also possible along some propagation directions.
To develop a highorder compact ADI FD scheme, we apply the Padé approximation to the secondorder central FD operators, so the second derivatives and can approximated with fourthorder accuracy. For example,
(8)  
(9) 
The fourthorder approximations of and can be obtained similarly.
Let , , . Substituting the fourthorder Padé approximations into Eq. (1) gives
(10) 
Truncation error analysis shows that the algorithm is fourthorder accurate in time and space with the truncation error , provided the solution and satisfy certain smoothness conditions. As shown in Liao2018 , the difficulty to develop highorder compact scheme for wave equation with nonconstant velocity is that, one can not multiply the operator
(11) 
to both sides of Eq. (10) to cancel the fractional operators , , and , which are difficult to implement. To overcome this problem, we develop an algebraic strategy which will be described in the next section.
3 Derivation of the compact highorder ADI method
Now we extend the novel algebraic manipulation introduced in Liao2018 to the threedimensional acoustic wave equation with variable velocity. We first demonstrate the difficulty in applying the Padé approximation to the finite difference operator for solving the acoustic wave equation with nonconstant velocity. Multiplying the operator in Eq. (11) to Eq. (10) yields
(12) 
We use the first term on the righthand side to illustrate the issue here. Since , and are commutative, we change the order of the three finite difference operators to obtain
(13) 
As discussed previously, when is nonconstant in terms of , the operator and are not commutative. Hence, . Therefore, the operator does not cancel the operator . In other words,
(14) 
To solve this problem, a novel factorization technique is used to preserve the compactness and fourthorder convergence of the numerical scheme. Multiplying to Eq. (10) yields
(15) 
Collecting the term , we have
(16) 
Factoring the lefthand side of Eq. (16) yields
(17) 
where the factorization error is given by
(18) 
Using Taylor series, one can verify that , provided that and satisfy some conditions on smoothness. The result regarding the order of the truncation error is included in the following theorem.
Theorem 3.1

A detailed proof of the theorem is given in Appendix A. ∎

If Eq. (16) is factorized in a different order of , and , for instance as
(19) then the factoring error is given by
(20) which has the same error estimation as that defined in Eq. (
62).
Ignoring the factoring error ERR in Eq. (17) leads to the following compact fourthorder FD method
(21) 
Using ADI method, Eq. (21) can be efficiently solved in three steps
(22)  
(23)  
(24) 
Apparently the three equations are difficult to implement due to the three operators , and . To overcome this problem, we apply the following strategy. Firstly, divide both sides of Eq. (22) by , then multiply , we have
(25) 
Eq. (25) is still hard to implement because of the terms and . Substituting with , with , respectively, we obtain
(26) 
Note the difference between Eq. (25) and Eq. (26) is , therefore, the method is still fourthorder in space. It is worth to mention that the righthand side of Eq. (26) includes larger stencil in both and directions. Furthermore, larger stencil needs values of outside the boundary when and . To overcome this problem, we use onesided approximation to approximate the values outside of the boundary. For example, is approximated by a linear combination of with fourthorder accuracy. This boundary treatment is not complicate in terms of implementation, since it only involves the values at time level , which is known. Similarly, dividing then multiplying to both sides of Eq. (23) lead to
(27)  
Finally, Eq. (24) can be transformed to the equivalent linear system
Comments
There are no comments yet.