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
where are the displacement vector at location and time t, is the density, and 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
Let , and 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:
where 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 for it has the second-order accuracy . And when and for it has the fourth-order accuracy . The general analytical expression of is given in Appendix. Similarly, we can define the operators and . 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 is . In order to compute and and analyze the stability, we define the following functional spaces
where 0 represents the integer grid and the half integer grid or or . The other functional spaces , , , , can be defined similarly. For saving space, we omit their definitions. Let the scalar inner product be defined in by
Other inner products such as , and so on have similar meaning. In the following, we compute and respectively.
We have conservation of the discrete energy, that is: . The stability of the scheme will be proven if the potential energy and the kinetic energy are positive. Since is obviously positive, we need to find out under what conditions is positive.
The problem can be reformulated as: , and with
we look for the corresponding conditions. Since
we can bound I as follows:
Obviously, . In the following we estimate , and respectively. Since
we have the following estimates for :
Thus we obtain
Then for , we have
Similarly, we have
Substituting (21), (22) and (24) into (19), we have
where . 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
If , then it reduces to
where , and , and 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, amplitude, and
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 defined by
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
here and 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 as follow:
which is a function of , where indicates or 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 and 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 and azimuth 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 1. The normalized phase error for different CFL number at the stability limit . The propagation angles and are fixed.
Figure 2. The normalized phase error for different CFL number at different propagation angle . The propagation angle is fixed.
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 and , where k is the wave number, and 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
 De Basabe, J.D. and Sen, M.K. (2010) Stability of the High-Order Finite Elements for Acoustic or Elastic Wave Propagation with High-Order Time Stepping. Geophysical Journal International, 181, 577-590.
 Bécache, E., Joly, P. and Tsogka, C. (2002) A New Family of Mixed Finite Elements for the Linear Elastodynamic Problem. SIAM Journal on Numerical Analysis, 39, 2109-2132.
 Zhang, W., Chung, E.T. and Wang, C. (2014) Stability for Imposing Absorbing Boundary Conditions in the Finite Element Simulation of Acoustic Wave Propagation. Journal of Computational Mathematics, 32, 1-20.
 Komatitsch, D., Martin, R., Tromp, J., Taylor, M.A. and Wingate, B.A. (2001) Wave Propagation in 2-D Elastic Media Using a Spectral Element Method with Triangles and Quadrangles. Journal of Computational Acoustics, 9, 703-718.
 Dumbser, M., Käser, M. and Toto, E.F. (2007) An Arbitrary High-Order Discontinuous Galerkin Method for Elastic Waves on Unstructured Meshes-V. Local Time Stepping and p-Adaptivity. Geophysical Journal International, 171, 695-717.
 Dumbser, M., Käser, M. and de la Puente, J. (2007) Arbitrary High-Order Finite Volume Schemes for Seismic Wave Propagation on Unstructured Meshes in 2D and 3D. Geophysical Journal International, 171, 665-694.
 Zhang, W., Zhuang, Y. and Chung, E.T. (2007) A New Spectral Finite Volume Method for Elastic Wave Modelling on Unstructured Meshes. Geophysical Journal International, 206, 292-307.
 Bayliss, A., Jordan, K.E., Lemesurier, B. and Turkel, E. (1986) A Fourth Accurate Finite Difference Scheme for the Computation of Elastic Waves. Bulletin of the Seismological Society of America, 76, 1115-1132.
 Zingg, D.W., Lomax, H. and Jurgens, H. (1996) High-Accuracy Finite-Difference Schemes for Linear Wave Propagation. SIAM Journal on Scientific Computing, 17, 328-346.
 Cohen, G. and Joly, P. (1996) Construction and Analysis of Fourth-Order Finite Difference Schemes for the Acoustic Wave Equation in Nonhomogeneous Media. SIAM Journal on Numerical Analysis, 33, 1266-1302.
 Nilsson, S., Petersson, N.A., Sjögreen, B. and Kreiss, H.-O. (2007) Stable Difference Approximations for the Elastic Wave Equation in Second Order Formulation. SIAM Journal on Numerical Analysis, 45, 1902-1936.
 Zhang, W. (2019) A New Family of Fourth-Order Locally One-Dimensional Schemes for the 3D Elastic Wave Equation. Journal of Computational and Applied Mathematics, 348, 246-260.