Numerical simulation of wave propagation has important applications in many scientific fields such as geophysics and seismic inversion. There are several types of numerical methods to solve the wave equations, for example, the finite difference method, the finite element method     , the spectral element method  , the discontinuous Galerkin method   and the finite volume method   . Each of the above numerical methods has its own advantages and disadvantages. In this paper, we consider the finite difference method.
The finite difference method is a very popular method because of high computational efficiency. In fact, it has been applied to wave simulation for several decades     . Since perfect numerical simulation depends on both stability and the order of accuracy, the high-order schemes and the corresponding stability are an important research topic of this field. In particularly, we may list a few here. In  , Cohen and Joly construct and analyses a family of fourth-order schemes for the acoustic wave equation in nonhomogeneous media. In  , Sei analysis the stability of high-order difference schemes for the 2D elastic wave equation in heterogeneous media. The stable difference approximation for the 3D elastic wave equation in the second-order formulation in heterogeneous media has been investigated in  . In  , a new family of locally one-dimensional schemes with fourth-order accuracy both in space and time for the 3D elastic wave equation is constructed and the stability is derived. The constructed new schemes in  only involve a three-point stencil in each spatial direction to achieve fourth-order accuracy. In this paper, based on the energy method, we study the stability analysis for the high-order staggered-grid schemes of the 3D elastic wave equation in heterogeneous media. To our knowledge, there is no work in this respect and our result is new.
The reminder of the paper is organized as follows. In Section 2, we present the governing equation and the high-order difference schemes in heterogeneous on staggered-grid grids. In Section 3, the stability analysis for the high-order difference schemes in heterogeneous is presented. In Section 4, the plane wave analysis in homogeneous case is investigated. In Section 5, we present numerical comparisons for 3D elastic wave simulation. Finally the conclusion and discussions are given in Section 6.
2. High-Order Spatial Discretization
We consider the following three-dimensional (3D) elastic wave equations in isotropic heterogeneous media
are the displacement vector at location
and time t,
is the density,
are the Lamé parameters,
is the external force.
Using the stress tensor, we can formulate the above system (1) as a first order in the following ways
be the spatial steps of
directions respectively. Now discretization of (2) with the second-order accuracy in space gives
For computing this, we need the values of u, v, and w at the grid
. One convenient way is to choose averaging the corresponding vales. For example,
However, such choices have no physical meaning. Another way is to compute
directly on staggered grids. In particular, we replace Equations (4)-(6) with Equations (7)-(9):
Obviously, the schemes (7)-(9) are the second-order accuracy in space. In order to construct high-order accuracy scheme in space, we first define the following functional spaces. Now we introduce the differentiation operator
on half integer grids with
order as follows:
is difference coefficients on the staggered grids, which can be calculated by a fast algorithm    by Matlab tool. Obviously, the approximation (10) or (11) has
order accuracy. For example, when
it has the second-order accuracy
. And when
it has the fourth-order accuracy
. The general analytical expression of
is given in Appendix. Similarly, we can define the operators
. Here, the subscript of the operator refers to the direction of differentiation.
Now, we can construct the semi-discrete schemes of system (1)
Applying the central difference approximation for time with the second-order accuracy, we obtain the full-discrete schemes of system (1), we obtain
where n denotes the time index and
the time step.
3. Stability Analysis
We now turn to the study of the numerical stability of the schemes (15)-(17). We are going to proceed by the energy method in analogy with continuous energy given by:
In a source-free infinite medium, the energy is conservative, i.e.,
. The discrete energy at time
. In order to compute
and analyze the stability, we define the following functional spaces
where 0 represents the integer grid
the half integer grid
. The other functional spaces
can be defined similarly. For saving space, we omit their definitions. Let the scalar inner product be defined in
Other inner products such as
and so on have similar meaning. In the following, we compute
Thus we obtain
, we have
Similarly, we have
Substituting (21), (22) and (24) into (19), we have
. Thus we obtain the sufficient stability condition for the numerical scheme (15)-(17). If the grid is uniform, i.e.
, then (26) gives
Therefore we summarize the conclusion above into the following theorem.
Theorem 1. A sufficient stability condition for the numerical schemes (15)-(17) is
, then it reduces to
are given by (20), (23) and (25) respectively.
4. Plane Wave Analysis
We turn to Fourier analysis  and we will derive the dispersion relation and by the von Neumann criterion we will get a necessary and sufficient stability condition. In homogeneous case for (12)-(14), the full-discrete schemes can be written as
We assume that
is a solution of Equation (29)-(31), where
is the angular frequency,
is the wave vector. Here
is the propagation angle and
the propagation azimuth. The two angles determine the movement direction of the plane wave in the 3D space.
Substituting the plane wave solution into Equations (29)-(31), we obtain the following relations:
By introducing the matrix B with elements
we can write the relations (32)-(34) as the following matrix form
The eigenvalues of B then express
as a function of
, which is the dispersion relation. There are three eigenvalues for matrix B. One eigenvalue is corresponding to the longitudinal or compressional wave, the double eigenvalues are corresponding to the transverse or shear wave. Thus we have the following two different relations
are the velocities of compressional and shear waves. Note that
is always larger than
. With the dispersion relations (36), we can apply the von Neumann stability criterion. A necessary stability is that the eigenvalues of B must be lower than 1. Thus we have
It is easy to verify that
Therefore we obtain the following theorem.
Theorem 2. In the homogeneous case, a sufficient and necessary stability condition for the numerical schemes (29)-(31) is given by
Proof Combining Equations (37) and (38), we obtain (39). Moreover, the matrix B is symmetric in homogeneous case. So the condition (39) is a sufficient and necessary condition. The proof is complete.
We now define the normalized phase error
which is a function of
which is related to different kinds of compressional wave and shear wave.
The stability condition
is defined by Courant-Friedrichs-Lewy (CFL) condition which bounds the interval for stability. We plot some dispersion curves based on Equation (40). Without loss of generality, we present dispersion curves for some special propagation angle and azimuth. Figure 1 is the normalized phase error for fixed
with different values of CFL condition and it shows that the phase error drops as increasing the order of accuracy. Figure 2 shows the normalized phase error for
and different values of
for different order or L. The figures for other propagation angle
are similar we omit them for saving space.
5. Numerical Computations
Wave simulation ignited by a point source is usually adopted in geophysical applications. For convenience, we simulate 3D elastic wave propagation in a homogeneous cubic model. The computational domain is
. The source is located in the center of the model and its time function is given by
which is loaded on the u component. The compressional velocity is 4000 m/s and the shear velocity 2500 m/s. The time step is
and the space step is
. Figure 3 shows the 3D snapshot of u component at propagation time 0.2 s. For brevity, we present some 2D slices of the 3D snapshots of u, v, and w components. The xz sections of 3D snapshots of u, v, and w components at propagation time 0.2 s are shown in Figures 4-6 respectively. We omit other sections for space. In our computations the scheme with fourth-order accuracy in space is applied. We remark that the comparisons between the numerical solution and the exact solution can be found in  . From Figures 4-6, we can clearly see the two types of waves, i.e. the compressional wave and the shear wave, which is consistent with the physical phenomenon.
Figure 3. The 3D snapshot of u component at propagation time 0.2 s.
Figure 4. The xz-section of 3D snapshot of u component at propagation time 0.2 s.
Figure 5. The xz-section of 3D snapshot of v component at propagation time 0.2 s.
Figure 6. The xz-section of 3D snapshot of w component at propagation time 0.2 s.
The staggered-grid difference method is a very important technique to solve wave equations numerically because of its high efficiency and the character of energy preservation. It has been well applied to seismic wave propagation for more than two decades. Based on the energy estimate method, we implement the stability analysis for the high-order staggered-grid schemes of the inhomogeneous 3D elastic wave equation. The stability result is controlled by the space varying parameters and the difference coefficients. The plane wave analysis in homogeneous media is completed and by the von Neumann criterion a necessary and sufficient stability condition is obtained. The analysis is helpful to design the computational parameters such as the time step and the space steps. Numerical computations are given to verify the effectiveness of the schemes. The key point of this paper is the theoretical analysis. In the future, we will consider more numerical computations for inhomogeneous media.
This work is supported by National Natural Science Foundation of China under grant numbers 11471328 and 51739007. It is also partially supported by National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences.
Appendix: Expression of the Coefficient
The calculation of difference coefficients on both regular and staggered-grid grids has been investigated by several authors     . In this appendix we present the analytical expression of the difference coefficient
in (10) or (11). In order to calculate the coefficient
, we consider an explicit staggered-grid difference expression for a function
where h is the step size, L is a positive number and
are the difference coefficients. Now consider
, where k is the wave number,
is a constant then we have
Then Taylor’s series expansion gives
Now equating the coefficient of
both sides we get
We can rewrite this equation in the following form
Now solving the above system, we get the following solutions