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 non-zero point source function is included in the equation. In particular, the high-order FD methods have attracted the interests of many researchers working on seismic modelling (see Chen2007 ; Cohen1996 ; Liu2009a ; Takeuchi2000 and references therein) due to the high-order 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 high-order FD schemes for the acoustic equations, and many accurate and efficient methods have been developed. Levander Levander1988 addressed the cost-effectiveness of solving real problems using high-order 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 time-space domain FD scheme with error O for 1-D, 2-D and 3-D 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 fourth-order 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 . High-order FD method is of particular importance for large-scale 3D acoustic wave equation, as it requires less grid pointsEtgen2007 .
These methods are accurate but are non-compact, which give rise to two issues: efficiency and difficulty in boundary condition treatment. For example, the conventional non-compact fourth-order FD scheme requires a five-point stencil in 1D to approximate , while in 2D, a 9-point stencil is needed in approximating . In 3D problem, a 13-point stencil is required to approximate the derivative
with fourth-order accuracy. To overcome these difficulties, a variety of compact higher-order FD schemes have been developed for hyperbolic, parabolic and elliptical partial differential equations. InChu1998 , the authors developed a family of fourth-order three-point combined difference schemes to approximate the first- and second-order spatial derivatives. In Kim2007 , the authors introduced a family of three-level implicit FD schemes which incorporate the locally one-dimensional method. For more recent compact higher-order difference methods, the readers are referred to Shukla2005 .
For three-dimensional 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 three-dimensional problem into a sequence of one-dimensional 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 fourth-order 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 high-order 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 low-dispersion 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 fourth-order compact ADI FD scheme was proposed to solve the two-dimensional acoustic wave equation with spatially variable velocity. Here we extended this method and the Padé approximation based high-order compact FD scheme in Das2014 to the 3D acoustic wave equation with non-constant velocity. The new method is compact and efficient, with fourth-order 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 second-order central FD schemes, compact high-order FD schemes for spatial derivatives and some other related high-order method in Section 2, then derive the new compact fourth-order 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
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 second-order central FD schemes are then given by approximating the second derivatives in Eq.(1). For example,
The approximations of and are similar to Eq. (6).
To improve the method to fourth-order in space, the conventional high-order 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
which can be as accurate as -order in . The conventional high-order 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 time-domain high-order FD methods have been derived by Liu and Sen Liu2009a . The idea of the time-domain high-order FD method is to determine coefficients using time-space domain dispersion. As a result, the coefficient will be a function of . It was noted that in 1D case, the time-domain high-order 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 high-order compact ADI FD scheme, we apply the Padé approximation to the second-order central FD operators, so the second derivatives and can approximated with fourth-order accuracy. For example,
The fourth-order approximations of and can be obtained similarly.
Let , , . Substituting the fourth-order Padé approximations into Eq. (1) gives
Truncation error analysis shows that the algorithm is fourth-order 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 high-order compact scheme for wave equation with non-constant velocity is that, one can not multiply the operator
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 high-order ADI method
Now we extend the novel algebraic manipulation introduced in Liao2018 to the three-dimensional 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 non-constant velocity. Multiplying the operator in Eq. (11) to Eq. (10) yields
We use the first term on the right-hand side to illustrate the issue here. Since , and are commutative, we change the order of the three finite difference operators to obtain
As discussed previously, when is non-constant in terms of , the operator and are not commutative. Hence, . Therefore, the operator does not cancel the operator . In other words,
To solve this problem, a novel factorization technique is used to preserve the compactness and fourth-order convergence of the numerical scheme. Multiplying to Eq. (10) yields
Collecting the term , we have
Factoring the left-hand side of Eq. (16) yields
where the factorization error is given by
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.
A detailed proof of the theorem is given in Appendix A. ∎
Ignoring the factoring error ERR in Eq. (17) leads to the following compact fourth-order FD method
Using ADI method, Eq. (21) can be efficiently solved in three steps
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
Eq. (25) is still hard to implement because of the terms and . Substituting with , with , respectively, we obtain
Note the difference between Eq. (25) and Eq. (26) is , therefore, the method is still fourth-order in space. It is worth to mention that the right-hand 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 one-sided approximation to approximate the values outside of the boundary. For example, is approximated by a linear combination of with fourth-order 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
Finally, Eq. (24) can be transformed to the equivalent linear system