doi:10.1155/2008/291968 Research Article Higher-Order Splitting Method for Elastic Wave Propagation

Motivated by seismological problems, we have studied a fourth-order split scheme for the elastic wave equation. We split in the spatial directions and obtain locally one-dimensional systems to be solved. We have analyzed the new scheme and obtained results showing consistency and stability. We have used the split scheme to solve problems in two and three dimensions. We have also looked at the influence of singular forcing terms on the convergence properties of the scheme.


Introduction
We are motivated by seismological problems that can be studied with higher-order splitting methods scheme.Because of the decomposition, we can save memory and computational time, which is important to study realistic elastic wave propagation.The ideas behind are to split in the spatial directions and obtain locally one-dimensional systems to be solved.Traditionally, using the classical operator splitting methods, we decouple the differential equation into more basic equations.These methods are often not sufficiently stable while also neglecting the physical correlations between the operators.Inspired by the work for the scalar wave equation presented in 1 , we devise a fourth-order split scheme for the elastic wave equation.From there on, we are going to develop new efficient methods based on a stable variant by coupling new operators and deriving new strong directions.We are going to examine the stability and consistency analyses for these methods and adopt them to linear acoustic wave equations seismic waves .Numerical experiments can validate our theoretical results and show the possibility to apply our methods.
The paper is organized as follows.A mathematical model based on the wave equation is introduced in Section 2. The utilized discretization methods are described in Section 3. The splitting method for the scalar and vectorial wave equations are discussed in Section 4 and the stability and consistency analyses are given.We discuss the numerical experiments in International Journal of Mathematics and Mathematical Sciences Section 5 with respect to scalar and vectorial problems.Finally, in Section 6, we foresee our future works in the area of splitting and decomposition methods.

Mathematical model
The mathematical models are studied in the following subsection.We introduce a scalar and also a vectorial model to distinguish the splitting methods.

Scalar wave equation
The motivation for the study presented below is coming from a computational simulation of earthquakes, see 2 , and the examination of seismic waves, see 3, 4 .We concentrate on the scalar wave equation, see 1 , for which the mathematical equations are given by ∂ tt u D∇•∇u, in Ω, u x, 0 u 0 x , u t x, 0 u 1 x , in Ω.

2.1
The unknown function u u x, t is considered to be in Ω × 0, T ⊂ R d × R, where the spatial dimension is given by d.
For three dimensions, we define the diffusion tensor as which describes the wave propagation.Further, the diffusion tensor D is given anisotropic, with D 1 , D 2 , D 3 ∈ R for D 1 , D 2 ≥ D 3 .The functions u 0 x and u 1 x are the initial conditions for the wave equation.We deal with the following boundary conditions: where all boundary conditions are on ∂Ω × T .

Elastic wave propagation
We consider the initial-value problem for the elastic wave equation for constant coefficients, given as where U ≡ U x, t is the displacement vector with components u, v T or u, v, w T in two and three dimensions, f, g 0 , and g 1 are known initial functions, and finally x x, y, z T .This equation is commonly used to simulate seismic events.
In seismology, it is common to use spatially singular forcing terms which can look like where F is a constant vector.A numeric method for 2.4a needs to approximate the Dirac function correctly in order to achieve full convergence.

Discretization methods
In this section, we discuss the discretization methods, both for time and space, to construct higher order methods.Because of the combination of both discretization, we can further show also higher-order methods for the splitting schemes, see also 1 .

Discretization of the scalar equation
At first, we underly finite difference schemes for the time and spatial discretization.
For the classical wave equation, this discretization is the well-known discretization in time and space.
Based on this discretization, the time is discretized as where the index i refers to the space point x i and Δt t n 1 − t n is the time step.We apply finite difference methods for the spatial discretization.The spatial terms and the initial conditions are given as where the index n refers to the time t n and Δx x i 1 − x i is the grid width.

International Journal of Mathematics and Mathematical Sciences
Then the two-dimensional equation, is discretized with the unconditionally stable implicit η-method, see 5 , where Δx and Δy are the grid width in x and y and 0 ≤ η ≤ 1.The initial conditions are given by U x, y, t n u x, y, t n and U x, y, t n−1 u x, y, t n − Δtu t x, y, t n .These discretization schemes are adopted to the operator splitting schemes.On the finite differences grid, Δt corresponds to the time step, and h x , h y are the grid sizes in the different spatial directions.The time nΔt is denoted by t n , and i, j refer to the spatial coordinates of the grid point ih x , jh y .Let u n denote the grid function on the time level n, and u n i,j be the specific value of u n at point i, j.In Section 3.2, we describe the traditional splitting methods for the wave equation.

Discretization of the vectorial equation
One of the first practical difference scheme with central differences used everywhere was introduced in 3 .To save space we exemplify it and some newer schemes in two dimensions first.
If we discretize uniformly in space and time on the unit square, we get a grid with grid points x j jh, y k kh, t n nΔt, where h > 0 is the spatial grid size and Δt the time step.Defining the grid function U n j,k U x j , y k , t n , the basic explicit scheme is ρ where M 2 is a difference operator and we use the standard difference operator notation: M 2 is a second-order difference approximation of the right-hand side operator of 2.4a .This explicit scheme is stable for time steps satisfying 6 Replacing M 2 with M 4 , a fourth-order difference operator given by and using the modified equation approach to eliminate the lower-order error terms in the time difference 6 , we obtain the explicit fourth-order scheme ρ where M 2 2 is a second-order approximation to the squared right-hand side operator in 2.4a .As it only needs to be second-order accurate, M 2 2 has the same extent in space as M 4 and no more grid points are used.This scheme has the same time step restriction as 3.8 .
In 1 the following implicit scheme for the scalar wave equation was introduced:

International Journal of Mathematics and Mathematical Sciences
When θ 1/12, the error of this scheme is fourth-order in time and space.For this θ value, it is, however, only conditionally stable, allowing a time step approximately 45% larger than 3.8 for θ ∈ 0.25, 0.5 , it is unconditionally stable .
In order to make it competitive with the explicit scheme 3.10 , we provide an operator split version of the implicit scheme 3.11 .This is made complicated by the presence of the mixed derivative terms that couple different coordinate directions.

Higher-order splitting method for the wave equations
In this section, we discuss the splitting methods for the wave equations.The higher order results as a combination between the spatial and time discretization method and the weighting factors in the splitting schemes.

Traditional splitting methods for the scalar wave equation
Our classical method is based on the splitting method of 5, 7 .
The classical splitting methods alternating direction methods ADIs are based on the idea of computing the different directions of the given operators.Each direction is computed independently by solving more basic equations.The result combines all the solutions of the elementary equations.So we obtain more efficiency by decoupling the operators.
The classical splitting method for the wave equation starts from where the initial functions u 0 and u 1 are given.We could also apply for u 1 that u t n u t n − u t n−1 /Δt O Δt u 1 .Consequently, we have u t n−1 ≈ u 0 − Δtu 1 .The right-hand side f t is given as a force term.
The spatial discretization terms are given by where the approximated discretization is given by

4.3
We could decouple the equation into 3 simpler equations obtaining a method of second order:

4.4c
where the result is given as u t n 1 with the initial conditions u t n u 0 , u t n−1 u 0 − Δtu 1 , and η ∈ 0, 0.5 .A fully coupled method is given for η 0 and for 0 < η ≤ 1 the decoupled method consists of a composition of explicit and implicit Euler methods.
We have to compute the first equation 4.4a and get the result u that is a further initial condition for the second equation 4.4b ; after whose computation we obtain u.In the third equation 4.4c , we have to put u as a further initial condition and get the result u t n 1 .
The underlying idea consists of the approximation of the pairwise operators: which we can raise to second order.

Boundary splitting method for the scalar wave equation
The time-dependent boundary conditions also have to be taken into account for the splitting method.Let us consider the three-operator example with the equations

4.6
where A D 1 x, y, z ∂ 2 /∂x 2 , B D 2 x, y, z ∂ 2 /∂y 2 , and C D 3 x, y, z ∂ 2 /∂z 2 are the spatial operators.The wave-propagation functions are as follows: Hence, for 3 operators, we have the following second-order splitting method: where the result is given as u t n 1 .
The boundary values are given by the following.
is split in where n is the outer normal vector and the anisotropic diffusion D, see 2.2 , is the parameter matrix to the wave-propagations.
We have the following initial conditions for the three equations:

4.17
International Journal of Mathematics and Mathematical Sciences Remark 4.1.By solving the two or three splitting steps, it is important to mention that each solution u, u, and u is corrected only once by using the boundary conditions.Otherwise, an "overdoing" of the boundary conditions takes place.

LOD method: locally one-dimensional method for the scalar wave equation
In the following, we introduce the LOD method as an improved splitting method while using prestepping techniques.
The method was discussed in 1 and is given by

4.18
where η ∈ 0.0, 0.5 and A, B are the spatial discretized operators.
If we eliminate the intermediate values in 4.18 , we obtain where B η η 2 Δt 4 AB and thus B η u n 1 − 2u n u n−1 O Δt 6 .So, we obtain a higher-order method.
Remark 4.2.For η ∈ 0.25, 0.5 , we have unconditionally stable methods and for higher order we use η 1/12.Then, for sufficiently small time steps, we get a conditionally stable splitting method.

Stability and consistency analysis for the LOD method of the scalar wave equation
The consistency of the fourth-order splitting method is given in the next theorem.Hence, we assume discretization orders of O h p , p 2, 4, for the discretization in space, where h h x h y is the spatial grid width.
Then we obtain the following consistency result for our method 4.18 .
Theorem 4.3.The consistency of the LOD method is given by where ∂ tt is a second-order discretization in time and A is the discretized fourth-order spatial operator.
Proof.We add 4.18 and obtain the following, see also 1 : Therefore, we obtain a splitting error of B u n 1 − 2u n u n−1 .Sufficient smoothness assumed that we have u n 1 − 2u n u n−1 O Δt 2 , and we obtain B u n 1 − 2u n u n−1 O Δt 4 .Thus, we obtain a fourth-order method if the spatial operators are also discretized as fourth-order terms.
The stability of the fourth-order splitting method is given in the following theorem.
Proof.We have to proof the theorem for a test function ∂ t u n , where ∂ t denotes the central difference.
For j ≥ 1, we have Multiplying with Δt and summarizing over j yield We can derive the identities

4.25
and obtain the result see also the idea of 1 .
To compute the error of the local splitting, we have to use the multiplier A 1 A 2 , thus for large constants, we have an unconditional small time step.Remark 4.6. 1 The unconditional stable version of LOD method is given for θ ∈ 0.25, 0.5 .
3 θ 1/12, we have a fourth-order method in time O Δt 2 h p , p ≥ 2. 4 θ 0 we have a second-order explicit scheme.5 For the stable version of the LOD method, the CFL condition should be taken into account for all θ ∈ 0, 0.25 with CFL Δt 2 /Δx 2 max D max , where x max are the maximal spatial step and D max are the maximal wave-propagation parameter in space.
In the next subsections, we discuss the higher-order splitting methods for the vectorial wave equations.

Higher-order splitting method for the vectorial wave equation
In the following, we present a fourth-order splitting method based on the basic scheme 3.11 .We split the operator M 4 into three parts: M xx , M yy , and M xy , where we have

4.27
Our proposed split method has the following steps:

4.28
Here, the first step is explicit, while the second and third steps treat the derivatives along the coordinate axes implicitly and the mixed derivatives explicitly.This is similar to how the mixed case is handled for parabolic problems 8 .
Note that each of the equation systems that needs to be solved in steps 2 and 3 is actually two decoupled tri-diagonal systems that can be solved independently.

Stability and consistency of the higher-order splitting method of the vectorial wave equations
The consistency of the fourth-order splitting method is given in the following theorem.
We have for all sufficiently smooth functions U x, t the following discretization order:

4.29
Furthermore, the split operators are also discretized with the same order of accuracy.Then, we obtain the following consistency result for the split method 4.28 .
Theorem 4.7.The split method has a splitting error which for smooth solutions U is O Δt 4 , where it is assumed that Δt O h .
Proof.We assume in the following that f 0, 0 T .We add 4.28 and obtain, like in the scalar case 1 , the following result for the discretized equations O Δt 2 and we obtain It is important that the influence of the mixed terms can be also be discretized as fourth-order method and, therefore, the terms are canceled in the proof.
For the stability, we have to denote an appropriate norm, which is in our case the L 2 Ω .
In the following, we introduce the notation of the norms.
Remark 4.8.For our system, we extend the L 2 -norm as where U 2 u 2 v 2 or U 2 u 2 v 2 w 2 in two and three dimensions.
Remark 4.9.For a discrete grid function U n j,k , the L 2 -norm is given as where Δx is the uniform grid length in x and y, M is the number of grid points in x and y.Further, U n j,k is the solution at grid point j, k and at time t n .

Numerical experiments
In this section, we present the numerical experiments for scalar and vectorial wave equations.The benefit of the splitting methods is discussed.

Numerical examples of the scalar wave equation
To test examples for the scalar wave equations, we discussed numerical experiments, which are based on analytical solutions.We present various boundary conditions and also spatial-dependent propagation functions.The benefit of the splitting method to reduce the computational amount is discussed with respect to the approximation errors.

Test example 1: problem with analytical solution and Dirichlet boundary condition
We deal with a two-dimensional example with constant coefficients where we can derive an analytical solution: where the initial conditions can be written as u x, y, t n u 0 x, y and u x, y, t n−1 u x, y, t n 1 u x, y, Δt .The analytical solution is given by For the approximation error, we choose the L 1 -norm.
The L 1 -norm is given by err L 1 : i,j 1,...,m V i,j u x i , y j , t n − u ana x i , y j , t n , 5.3 where u x i , y j , t n is the numerical and u ana x i , y j , t n is the analytical solution and V i,j ΔxΔy.
Our test examples are organized as follows.
1 The non-stiff case.We choose D 1 D 2 1 with a rectangle as our model domain Ω 0, 1 × 0, 1 .We discretize with Δx 1/16 and Δy 1/16 and Δt 1/32 and choose our parameter η between 0 ≤ η ≤ 1.The exemplary function values u num and u ana are taken from the center of our domain.2 The stiff case.We choose D 1 D 2 0.01 with a rectangle as our model domain Ω 0, 1 × 0, 1 .We discretize with Δx 1/32 and Δy 1/32 and Δt 1/64 and choose our parameter η between 0 ≤ η ≤ 1.The exemplary function values u num and u ana are taken from the point 0.5, 0.5625 .
The experiments are done with the uncoupled standard discretization method i.e., the finite differences methods for time and space, and with the operator splitting methods, i.e., the classical operator splitting method and the LOD method .
The numerical errors for the non-stiff case with Dirichlet boundary conditions are presented in Figure 1 and their results in Figure 2. The numerical errors for the stiff case with Dirichlet boundary conditions are presented in Figure 3 and their results in Figure 4.
Remark 5.1.In the experiments, we compare the non-splitting with the splitting methods.We obtain nearly the same results and could see improved results for the LOD method, which is for η 1/12 a fourth-order method.

Test example 2: problem with analytical solution and Neumann boundary condition
In this example, we modify our boundary conditions with respect to the Neumann boundary.We deal with our 2D example where we can derive an analytical solution: where Ω 0, 1 × 0, 1 .D 1 1, D 2 0.5 and the initial conditions can be written as u x, y, t n u 0 x, y and u x, y, t n−1 u x, y, t n 1 u x, y, Δt .The analytical solution is given as We have the same discretization methods as in test example 1.
The numerical errors for the non-stiff case with Neumann boundary conditions are presented in Figure 5 and their results in Figure 6.Remark 5.2.In the experiments, we can obtain the same accuracy as for the Dirichlet boundary conditions.More accurate results are gained by the LOD method with small η.We obtain also stable results in our computations.

Test example 3: spatial-dependent wave equation
In this experiment, we apply our method to the spatial-dependent problem, given by where D 1 x, y 0.1x 0.01y 0.01, D 2 x, y 0.01x 0.1y 0.1.To compare the numerical results, we cannot use an analytical solution, that is why in a first prestep we are computing a reference solution.The reference solution is done with the finite difference scheme with fine time and space steps.
Concerning the choice of the time steps, it is important to consider the CFL condition, that is now based on the spatial coefficients.
The numerical errors for the spatial-dependent parameters with Dirichlet boundary conditions are presented in Figure 7 and their results in Figure 8.
The numerical errors for the spatial-dependent parameters with Neumann boundary conditions are presented in Figure 9 and their results in Figure 10.

Remark 5.4.
In the experiments, we analyze the classical operator splitting and the LOD method and show that the LOD method yields yet more accurate values.

Numerical experiments of the elastic wave propagation
To test a fourth-order split method, we have done grid convergence studies on two types of problems.For the first, we impose a smooth solution of 2.4a using a specific form of the forcing function f and check the error of the numerical solution against the known solution as the grid is refined.
For the second problem, we use a singular forcing function 2.5 , and compare the numerical solution to a solution computed using the Green's function for the free space elastodynamic problem.The convergence for this case is dependent not only on the approximations of time and space derivatives but also on how the Dirac function is approximated.

International Journal of Mathematics and Mathematical Sciences
During the numerical testing, we have observed a need to reduce the allowable time step when the ration of λ over μ became too large.This is likely from the influence of the explicitly treated mixed derivative.For really high ratios >20 , a reduction of 35% was necessary to avoid numerical instabilities.

Initial values and boundary conditions
In order to start the time stepping scheme, we need to know the values at two earlier time levels.Starting at time t 0, we know the value at level n 0 as U 0 g 0 .The value at level n −1 can be obtained by Taylor expansion as where we use and also for 5.9c and 5.9d ,

5.9f
We are not considering the boundary value problem in this paper and so will not be concerned with constructing proper difference stencils at grid points close to the boundaries of the computational domain.We have simply added a two-point-thick layer of extra-grid points at the boundaries of the domain and assigned the correct analytical solution at all points in the layer every time step.
Remark 5.5.For the Dirichlet boundary conditions, the splitting method see 4.28 conserves also the conditions.We can use for the 3 equations see 4.28 , so for U * , U * * , and for U n 1 , the same conditions.For the Neumann boundary conditions and other boundary conditions of higher order, we have also to split the boundary conditions with respect to the split operators, see 12 .

5.11
Using the split method we solved 2.4a on a domain x × y −11 × −11 up to t 2. We used two sets of material parameters; for the first case ρ, λ, and μ were all equal to 1, for the second case ρ and μ were 1 and λ was set to 14. Solving on four different grids with a refinement factor of two in each direction between the successive grids we obtained the results shown in Table 1.The errors are measured in the ∞-norm defined as U j,k max max j,k |u j,k |, max j,k |v j,k | .As can be seen we get the expected 4th order convergence for problems with smooth solutions.
To check the influence of the splitting error N 4,θ on the error we solved the same problems using the non-split scheme 3.11 .The results are shown in Table 2.The errors are only marginally smaller than for the split scheme.

Singular forcing terms
In seismology and acoustics it is common to use spatially singular forcing terms which can look like International Journal of Mathematics and Mathematical Sciences f Fδ x g t , 5.12 where F is a constant direction vector.A numeric method for 2.4a needs to approximate the Dirac function correctly in order to achieve full convergence.Obviously we cannot expect convergence close to the source as the solution will be singular for two and three dimensional domains.The analyzes in 13, 14 demonstrate that it is possible to derive regularized approximations of the Dirac function, which result in point wise convergence of the solution away from the sources.Based on these analyzes, we define one 2nd δ h 2 and one 4th δ h 4 order regularized approximations of the one dimensional Dirac function, where in the above x x/h.The two and three dimensional Dirac functions are then approximated as δ h 2,4 x δ h 2,4 y and δ h 2,4 x δ h 2,4 y δ h 2,4 z .The chosen time dependence was a smooth function given by g t which is C ∞ .Using this forcing function we can compute the analytical solution by integrating the Green's function given in 15 in time.The integration was done using numerical quadrature routines from Matlab.Figures 11 and 12 shows examples of what the errors look like on a radius passing through the singular source at time t 0.8 for different grid sizes h and the two approximations δ h 2 and δ h 4 .As can be seen the error is smooth and converges a small distance away from the source.However, using δ h 2 limits the convergence to 2nd order, while using δ h 4 gives the full 4th order convergence away from the singular source.When t > 1 the forcing goes to zero and the solution will be smooth everywhere.Table 3 shows the convergence behavior at time t 1.1 for four different grids.Note that the full convergence is achieved even if the lower order δ h 2 is used as an approximation for the smooth material properties, making the two methods roughly comparable in computational cost.

Three-dimensional test example for the elastic wave propagation
The motivation to compute also the three dimensional elastic wave propagation arose from the need to understand the anisotropy of the different dimensions, see 2 .We apply the three-dimensional model 2.4a -2.4c for our proposed splitting schemes.

The splitting scheme
In three dimensions a 4th order difference approximation of the right hand side operator becomes operating on grid functions U n j,k,l defined at grid points x j , y k , z l , t n similarly to the two dimensional case.We can split M 4 into six parts; M xx , M yy , M zz containing the three second order directional difference operators, and M xy , M yz , M xz containing the mixed difference operators.method.We have designed with higher order spatial and time discretization methods the stabile higher order splitting methods.The benefit of the splitting methods is due to the different scales and therefore the computational process in decoupling the stiff and the nonstiff operators into different equation is accelerated.The LOD method as a 4th-oder method has the advantage of higher accuracy and can be used for such decoupling regards.For our realistic application in elastic wave propagation, the split scheme has been proven to work well in practice for different types of material properties It is comparable to the fully explicit 4th order scheme 3.10 in terms of computational cost, but should be easier to implement, as no difference approximations of higher order operators are needed.
In a next work we will discuss a model in seismology, which have to be more accurate in the boundary conditions.For such models we have to develop higher order stable splitting methods.

Figure 1 : 5 bFigure 2 :
Figure 1: Numerical error for standard and modified methods, with respect to the η parameter and Dirichlet boundary conditions.

Figure 5 :
Figure 5: Numerical error for standard and modified methods, with respect to the η parameter and Neumann boundary conditions.

InternationalFigure 7 : 1 bFigure 8 :
Figure 7: Numerical error for standard and modified methods, with respect to the η parameter, spatialdependent parameters, and Dirichlet boundary conditions.

Figure 9 : 5 bFigure 10 :
Figure 9: Numerical error for standard and modified methods with respect to the η parameter, spatialdependent parameters, and Neumann boundary conditions.
In the next test example, we study the Neumann boundary conditions.

Table 1 :
Errors in max-norm for decreasing h and smooth analytical solution U true .Convergence rate indicates a fourth-order convergence for the split scheme.

Table 2 :
Errors in max-norm for decreasing h and smooth analytical solution U true and using the non-split scheme.Comparing with Table1, we see that the splitting error is very small for this case.
f sin t − x sin y − 2μ sin t − x sin y − λ μ cos x cos t − y sin t − x sin y ,sin t − y sin x − 2V s 2 sin x sin t − y − λ μ cos t − x cos y sin y sin t − y T , true sin x − t sin y , sin y − t sin x T .

Table 3 :
Errors in max-norm for decreasing h and analytical solution U true .Convergence rate approaches 4th order after the singular forcing term goes to zero.

Table 4 :
Errors in max-norm for decreasing h and smooth analytical solution U true .Convergence rate indicates 4th order convergence for the three dimensional split scheme.