A compact high order Alternating Direction Implicit method for three-dimensional acoustic wave equation with variable coefficient

by   Keran Liao, et al.
University of Calgary

Efficient and accurate numerical simulation of seismic wave propagation is important in various Geophysical applications such as seismic full waveform inversion (FWI) problem. However, due to the large size of the physical domain and requirement on low numerical dispersion, many existing numerical methods are inefficient for numerical modelling of seismic wave propagation in a heterogeneous media. Despite the great efforts that have been devoted during the past decades, it still remains a challenging task in the development of efficient and accurate finite difference method for the multi-dimensional acoustic wave equation with variable velocity. In this paper, we proposed a Padé approximation based finite difference scheme for solving the acoustic wave equation in three-dimensional heterogeneous media. The new method is obtained by combining the Padé approximation and a novel algebraic manipulation. The efficiency of the new algorithm is further improved through the Alternative Directional Implicit (ADI) method. The stability of the new algorithm has been theoretically proved by the energy method. The new method is conditionally stable with a better Courant - Friedrichs - Lewy condition (CFL) condition, which has been verified numerically. Extensive numerical examples have been solved, which demonstrated that the new method is accurate, efficient and stable.



There are no comments yet.


page 20

page 21

page 22


An Efficient and high accuracy finite-difference scheme for the acoustic wave equation in 3D heterogeneous media

Efficient and accurate numerical simulation of 3D acoustic wave propagat...

Efficient and Stable Finite Difference Modelling of Acoustic Wave Propagation in Variable-density Media

In this paper, we consider the development and analysis of a new explici...

Dirac Assisted Tree Method for 1D Heterogeneous Helmholtz Equations with Arbitrary Variable Wave Numbers

In this paper we introduce a method called Dirac Assisted Tree (DAT), wh...

Random noise attenuation on finite-difference wave propagation using fuzzy transform

Fuzzy Transform (F-transform) has been introduced as an approximation me...

A re-formulization of the transfer matrix method for calculating wave-functions in higher dimensional disordered open systems

We present a numerically stable re-formulization of the transfer matrix ...

A parallel hybrid implementation of the 2D acoustic wave equation

In this paper, we propose a hybrid parallel programming approach for a n...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 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. In

Chu1998 , 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.

Theorem 3.1

Assume that is the solution of the acoustic wave equation defined by Eqs. ( 1 - 4), and the coefficient satisfies the smoothness condition . Then the truncation error given in Eq. (18) satisfies

where and are the step sizes in time, and , respectively.

  • 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


    then the factoring error is given by


    which has the same error estimation as that defined in Eq. (


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