Higher-Order Splitting Method for Elastic Wave Propagation

Abstract

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.

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

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

2.1. 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 𝜕_𝑡𝑡𝑢=𝐃∇⋅∇𝑢,inΩ,𝑢(𝐱,0)=𝑢₀(𝐱),𝑢_𝑡(𝐱,0)=𝑢₁(𝐱),inΩ.(2.1) The unknown function 𝑢=𝑢(𝐱,𝑡) is considered to be in Ω×(0,𝑇)⊂ℝ^𝑑×ℝ, where the spatial dimension is given by 𝑑.

For three dimensions, we define the diffusion tensor as 𝐷𝐃=¹000𝐷₂000𝐷₃∈ℝ^3,+×ℝ^3,+,(2.2) which describes the wave propagation. Further, the diffusion tensor 𝐃 is given anisotropic, with 𝐷₁,𝐷₂,𝐷₃∈ℝ⁺ for 𝐷₁,𝐷₂≥𝐷₃. The functions 𝑢₀(𝐱) and 𝑢₁(𝐱) are the initial conditions for the wave equation.

We deal with the following boundary conditions: 𝑢(𝐱,𝑡)=𝑢₃,Dirichletboundarycondition,𝜕𝑢(𝐱,𝑡)𝜕𝑛=0,Neumannboundarycondition,𝐃∇𝑢(𝐱,𝑡)=𝑢_out,Outflowboundarycondition,(2.3) where all boundary conditions are on 𝜕Ω×𝑇.

2.2. Elastic Wave Propagation

We consider the initial-value problem for the elastic wave equation for constant coefficients, given as 𝜌𝜕_𝑡𝑡𝐔=𝜇∇²𝑈+(𝜆+𝜇)∇(∇⋅𝑈)+𝐟,(2.4a)𝐔(𝑡=0,𝐱)=𝐠₀𝜕(𝐱),(2.4b)_𝑡𝐔(𝑡=0,𝐱)=𝐠₁(𝐱),(2.4c) where 𝐔≡𝐔(𝐱,𝑡) is the displacement vector with components (𝑢,𝑣)^𝑇 or (𝑢,𝑣,𝑤)^𝑇 in two and three dimensions, 𝐟, 𝐠₀, and 𝐠₁ are known initial functions, and finally 𝐱=(𝑥,𝑦,𝑧)^𝑇. This equation is commonly used to simulate seismic events.

In seismology, it is common to use spatially singular forcing terms which can look like 𝐟=𝐅𝛿(𝐱)𝑔(𝑡),(2.5) where 𝐅 is a constant vector. A numeric method for (2.4a) needs to approximate the Dirac function correctly in order to achieve full convergence.

3. 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].

3.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 𝑈_{𝑡𝑡,𝑖}=𝑈_𝑖^𝑛+1−2𝑈^𝑛_𝑖+𝑈_𝑖^𝑛−1Δ𝑡²,𝑈(𝑡^𝑛)=𝑢(𝐱,𝑡),𝑈_𝑡(𝑡^𝑛)=𝑢_𝑡(𝐱,𝑡),(3.1) where the index 𝑖 refers to the space point 𝑥_𝑖 and Δ𝑡=𝑡^𝑛+1−𝑡^𝑛 is the time step. We apply finite difference methods for the spatial discretization. The spatial terms and the initial conditions are given as 𝑈_{𝑥𝑥,𝑛}=𝑈^𝑛_𝑖+1−2𝑈^𝑛_𝑖+𝑈^𝑛_𝑖−1Δ𝑥²,𝑈(𝑡^𝑛)=𝑢(𝐱,𝑡),𝑈_𝑡(𝑡^𝑛)=𝑢_𝑡(𝐱,𝑡),(3.2) where the index 𝑛 refers to the time 𝑡_𝑛 and Δ𝑥=𝑥_𝑖+1−𝑥_𝑖 is the grid width.

Then the two-dimensional equation, 𝑢_𝑡𝑡=𝐷₁𝑢_𝑥𝑥+𝐷₂𝑢_𝑦𝑦inΩ,𝑢(𝑥,𝑦,0)=𝑢₀(𝑥,𝑦),𝑢_𝑡(𝑥,𝑦,0)=𝑢₁(𝑥,𝑦),𝑢(𝑥,𝑦,𝑡)=𝑢₂on𝜕Ω,(3.3) is discretized with the unconditionally stable implicit 𝜂-method, see [5], 𝑈^𝑛+1_𝑖,𝑗−2𝑈^𝑛_𝑖,𝑗+𝑈^𝑛−1_𝑖,𝑗Δ𝑡²=𝐷₁Δ𝑥²(𝜂(𝑈^𝑛+1_𝑖+1,𝑗−2𝑈^𝑛+1_𝑖,𝑗+𝑈^𝑛+1_{𝑖−1,𝑗})+(1−2𝜂)(𝑈^𝑛_𝑖+1,𝑗−2𝑈^𝑛_𝑖,𝑗+𝑈^𝑛_{𝑖−1,𝑗})+𝜂(𝑈^𝑛−1_𝑖+1,𝑗−2𝑈^𝑛−1_𝑖,𝑗+𝑈^𝑛−1_{𝑖−1,𝑗}+𝐷))₂Δ𝑦²(𝜂(𝑈^𝑛_𝑖,𝑗+1−2𝑈^𝑛_𝑖,𝑗+𝑈^𝑛_{𝑖,𝑗−1})+(1−2𝜂)(𝑈^𝑛_𝑖,𝑗+1−2𝑈^𝑛_𝑖,𝑗+𝑈^𝑛_{𝑖,𝑗−1})+𝜂(𝑈^𝑛−1_𝑖,𝑗+1−2𝑈^𝑛−1_𝑖,𝑗+𝑈^𝑛−1_{𝑖,𝑗−1})),(3.4)where Δ𝑥 and Δ𝑦 are the grid width in 𝑥 and 𝑦 and 0≤𝜂≤1. The initial conditions are given by 𝑈(𝑥,𝑦,𝑡^𝑛)=𝑢(𝑥,𝑦,𝑡^𝑛) and 𝑈(𝑥,𝑦,𝑡^𝑛−1)=𝑢(𝑥,𝑦,𝑡^𝑛)−Δ𝑡𝑢_𝑡(𝑥,𝑦,𝑡^𝑛).

These discretization schemes are adopted to the operator splitting schemes.

On the finite differences grid, Δ𝑡 corresponds to the time step, and ℎ_𝑥, ℎ_𝑦 are the grid sizes in the different spatial directions. The time 𝑛Δ𝑡 is denoted by 𝑡^𝑛, and 𝑖,𝑗 refer to the spatial coordinates of the grid point (𝑖ℎ_𝑥,𝑗ℎ_𝑦). Let 𝐮^𝑛 denote the grid function on the time level 𝑛, and 𝑢^𝑛_𝑖,𝑗 be the specific value of 𝐮^𝑛 at point 𝑖,𝑗.

In Section 3.2, we describe the traditional splitting methods for the wave equation.

3.2. 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 𝑥_𝑗=𝑗ℎ, 𝑦_𝑘=𝑘ℎ, 𝑡_𝑛=𝑛Δ𝑡, where ℎ>0 is the spatial grid size and Δ𝑡 the time step. Defining the grid function 𝐔^𝑛_𝑗,𝑘=𝐔(𝑥_𝑗,𝑦_𝑘,𝑡_𝑛), the basic explicit scheme is𝜌𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘Δ𝑡²=ℳ₂𝐔^𝑛_𝑗,𝑘+𝐟^𝑛_𝑗,𝑘,(3.5)where ℳ₂ is a difference operatorℳ₂=(𝜆+2𝜇)𝐷^𝑥2+𝜇𝐷^𝑦2(𝜆+𝜇)𝐷^𝑥₀𝐷^𝑦₀(𝜆+𝜇)𝐷^𝑥₀𝐷^𝑦₀(𝜆+2𝜇)𝐷^𝑦2+𝜇𝐷^𝑥2,(3.6) and we use the standard difference operator notation: 𝐷^𝑥₊𝑣_𝑗,𝑘=1ℎ(𝑣_𝑗+1,𝑘−𝑣_𝑗,𝑘),𝐷^𝑥₋𝑣_𝑗,𝑘=𝐷^𝑥₊𝑣_{𝑗−1,𝑘},𝐷^𝑥₀=12(𝐷^𝑥₊+𝐷^𝑥₋),𝐷^𝑥2=𝐷^𝑥₊𝐷^𝑥₋(3.7)ℳ₂ 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]ℎΔ𝑡<𝜆+3𝜇.(3.8)

Replacing ℳ₂ with ℳ₄, a fourth-order difference operator given byℳ₄=ℎ(𝜆+2𝜇)(1−²𝐷12^𝑥2)𝐷^𝑥2ℎ+𝜇(1−²𝐷12^𝑦2)𝐷^𝑦2ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(𝜆+2𝜇)(1−²𝐷12^𝑦2)𝐷^𝑦2ℎ+𝜇(1−²𝐷12^𝑥2)𝐷^𝑥2,(3.9) 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 𝜌𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘Δ𝑡²=ℳ₄𝐔^𝑛_𝑗,𝑘+𝐟^𝑛_𝑗,𝑘+Δ𝑡²12(ℳ²₂𝐔^𝑛_𝑗,𝑘+ℳ₂𝐟^𝑛_𝑖,𝑗+𝜕_𝑡𝑡𝐟^𝑛_𝑖,𝑗),(3.10)where ℳ²₂ is a second-order approximation to the squared right-hand side operator in (2.4a). As it only needs to be second-order accurate, ℳ²₂ has the same extent in space as ℳ₄ 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:𝜌𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘Δ𝑡²=ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘)+𝜃𝐟^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐟^𝑛_𝑗,𝑘+𝜃𝐟^𝑛−1_𝑗,𝑘.(3.11) 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.

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

4.1. 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 𝜕_𝑡𝑡𝑢(𝑡)=(𝐴+𝐵+𝐶)𝑢(𝑡)+𝑓(𝑡),𝑡∈(𝑡^𝑛,𝑡^𝑛+1),𝑢(𝑡^𝑛)=𝑢₀,𝑢^′(𝑡^𝑛)=𝑢₁,(4.1) where the initial functions 𝑢₀ and 𝑢₁ are given. We could also apply for 𝑢₁ that 𝑢^′(𝑡^𝑛)=(𝑢(𝑡^𝑛)−𝑢(𝑡^𝑛−1))/Δ𝑡+𝑂(Δ𝑡)=𝑢₁. Consequently, we have 𝑢(𝑡^𝑛−1)≈𝑢₀−Δ𝑡𝑢₁. The right-hand side 𝑓(𝑡) is given as a force term.

The spatial discretization terms are given by 𝜕𝐴=²𝜕𝑥²𝜕,𝐵=²𝜕𝑦²𝜕,𝐶=²𝜕𝑧²,(4.2) where the approximated discretization is given by𝐴𝑢(𝑥,𝑦,𝑧)≈𝑢(𝑥+Δ𝑥,𝑦,𝑧)−2𝑢(𝑥,𝑦,𝑧)+𝑢(𝑥−Δ𝑥,𝑦,𝑧)Δ𝑥²,𝐵𝑢(𝑥,𝑦,𝑧)≈𝑢(𝑥,𝑦+Δ𝑦,𝑧)−2𝑢(𝑥,𝑦,𝑧)+𝑢(𝑥,𝑦−Δ𝑦,𝑧)Δ𝑦²,𝐶𝑢(𝑥,𝑦,𝑧)≈𝑢(𝑥,𝑦,𝑧+Δ𝑧)−2𝑢(𝑥,𝑦,𝑧)+𝑢(𝑥,𝑦,𝑧−Δ𝑧)Δ𝑧².(4.3)

We could decouple the equation into 3 simpler equations obtaining a method of second order: ^0𝑥002𝑑𝑢^{𝑐0𝑥002𝑑𝑐}−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐵𝑢(𝑡^𝑛)+Δ𝑡²𝐶𝑢(𝑡^𝑛)+Δ𝑡²(𝜂𝑓(𝑡^𝑛+1)+(1−2𝜂)𝑓(𝑡^𝑛)+𝜂𝑓(𝑡^𝑛−1)),(4.4a)_{0𝑥002𝑑𝑐}𝑢−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐵(𝜂𝑢^{0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐶𝑢(𝑡^𝑛)+Δ𝑡²(𝜂𝑓(𝑡^𝑛+1)+(1−2𝜂)𝑓(𝑡^𝑛)+𝜂𝑓(𝑡^𝑛−1)),(4.4b)𝑢(𝑡^𝑛+1)−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐵(𝜂𝑢^{0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐶(𝜂𝑢(𝑡^𝑛+1)+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²(𝜂𝑓(𝑡^𝑛+1)+(1−2𝜂)𝑓(𝑡^𝑛)+𝜂𝑓(𝑡^𝑛−1)),(4.4c) where the result is given as 𝑢(𝑡^𝑛+1) with the initial conditions 𝑢(𝑡^𝑛)=𝑢₀, 𝑢(𝑡^𝑛−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 ^0𝑥002𝑑𝑢^{𝑐0𝑥002𝑑𝑐} that is a further initial condition for the second equation (4.4b); after whose computation we obtain ^0𝑥0𝑢^02𝑑𝑐. In the third equation (4.4c), we have to put ^0𝑥0𝑢^02𝑑𝑐 as a further initial condition and get the result 𝑢(𝑡^𝑛+1).

The underlying idea consists of the approximation of the pairwise operators: Δ𝑡²𝐴𝜂(𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1))≈0,Δ𝑡²𝐵𝜂(𝑢^{0𝑥002𝑑𝑐}−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1))≈0,(4.5) which we can raise to second order.

4.2. 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 𝜕_𝑡𝑡𝑢(𝑡)=(𝐴+𝐵+𝐶)𝑢(𝑡)+ℎ(𝑡),𝑡∈(𝑡^𝑛,𝑡^𝑛+1),𝑢(𝑡^𝑛)=𝑔(𝑡),𝑢^′(𝑡^𝑛)=𝑓(𝑡),(4.6) where 𝐴=𝐷₁(𝑥,𝑦,𝑧)(𝜕²/𝜕𝑥²), 𝐵=𝐷₂(𝑥,𝑦,𝑧)(𝜕²/𝜕𝑦²), and 𝐶=𝐷₃(𝑥,𝑦,𝑧)(𝜕²/𝜕𝑧²) are the spatial operators. The wave-propagation functions are as follows: 𝐷₁(𝑥,𝑦,𝑧),𝐷₂(𝑥,𝑦,𝑧),𝐷₃(𝑥,𝑦,𝑧)∶ℝ³→ℝ⁺.(4.7)

Hence, for 3 operators, we have the following second-order splitting method: ^0𝑥002𝑑𝑢^{𝑐0𝑥002𝑑𝑐}−2𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛)+𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛)+𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛−1))+Δ𝑡²𝐵𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛)+Δ𝑡²𝐶𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}(𝑡^𝑛)+Δ𝑡²(𝜂ℎ(𝑡^𝑛+1)+(1−2𝜂)ℎ(𝑡^𝑛)+𝜂ℎ(𝑡^𝑛−1)),_{0𝑥002𝑑𝑐}𝑢−2𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛)+𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛)+𝜂𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛−1))+Δ𝑡²𝐵(𝜂𝑢^{0𝑥002𝑑𝑐}+(1−2𝜂)𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛)+𝜂𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛−1))+Δ𝑡²𝐶𝑢^{0𝑥002𝑑𝑐}(𝑡^𝑛)+Δ𝑡²(𝜂ℎ(𝑡^𝑛+1)+(1−2𝜂)ℎ(𝑡^𝑛)+𝜂ℎ(𝑡^𝑛−1)),𝑢(𝑡^𝑛+1)−2𝑢(𝑡^𝑛)+𝑢(𝑡^𝑛−1)=Δ𝑡²𝐴(𝜂𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐵(𝜂𝑢^{0𝑥002𝑑𝑐}+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²𝐶(𝜂𝑢(𝑡^𝑛+1)+(1−2𝜂)𝑢(𝑡^𝑛)+𝜂𝑢(𝑡^𝑛−1))+Δ𝑡²(𝜂ℎ(𝑡^𝑛+1)+(1−2𝜂)ℎ(𝑡^𝑛)+𝜂ℎ(𝑡^𝑛−1)),(4.8) where the result is given as 𝑢(𝑡^𝑛+1).

The boundary values are given by the following.

(1) Dirichlet values. We have to use the same boundary values for all 3 equations.

(2) Neumann values. We have to decouple the values into the different directions: 𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝜕𝑛=0(4.9) is split in𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛𝜕𝑥^𝑥+𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛𝜕𝑦^𝑦+𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛𝜕𝑧^𝑧𝜕=0;_{0𝑥002𝑑𝑐}𝑢𝜕𝑛=0(4.10) is split in 𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛𝜕𝑥^𝑥+𝜕𝑢^{0𝑥002𝑑𝑐}𝑛𝜕𝑦^𝑦+𝜕𝑢^{0𝑥002𝑑𝑐}𝑛𝜕𝑧^𝑧=0;𝜕𝑢(𝑡^𝑛+1)𝜕𝑛=0(4.11) is split in 𝜕𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛𝜕𝑥^𝑥+𝜕𝑢^{0𝑥002𝑑𝑐}𝑛𝜕𝑦^𝑦+𝜕𝑢^𝑛+1𝑛𝜕𝑧^𝑧=0.(4.12)

(3) Outflowing values, we have to decouple the values into the different directions: 𝐧𝐃∇𝑢^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}=𝑢_out(4.13) is split in 𝐷₁𝜕𝑢_𝑥^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛_𝑥+𝐷₂𝜕𝑢_𝑦^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛_𝑦+𝐷₃𝜕𝑢_𝑧^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛_𝑧=𝑢_out;𝐧𝐃∇𝑢^{0𝑥002𝑑𝑐}=𝑢_out(4.14) is split in 𝐷₁𝜕𝑢_𝑥^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛_𝑥+𝐷₂𝜕𝑢_𝑦^{0𝑥002𝑑𝑐}𝑛_𝑦+𝐷₃𝜕𝑢_𝑧^{0𝑥002𝑑𝑐}𝑛_𝑧=𝑢_out;𝐧𝐃∇𝑢^𝑛+1=𝑢_out(4.15) is split in 𝐷₁𝜕𝑢_𝑥^{0𝑥002𝑑𝑐0𝑥002𝑑𝑐}𝑛_𝑥+𝐷₂𝜕𝑢_𝑦^{0𝑥002𝑑𝑐}𝑛_𝑦+𝐷₃𝜕_𝑧𝑢^𝑛+1𝑛_𝑧=𝑢_out,(4.16) where 𝐧 is the outer normal vector and the anisotropic diffusion 𝐃, see (2.2), is the parameter matrix to the wave-propagations.

We have the following initial conditions for the three equations: 𝑢(𝑡^𝑛)=𝑢₀,𝑢(𝑡^𝑛−1)=𝑢₀−Δ𝑡𝑢₁+Δ𝑡²2((𝐴+𝐵+𝐶)𝑢₀+𝑓(𝑡^𝑛))+𝑂(Δ𝑡³),𝑢(𝑡^𝑛−1)=𝑢₀−Δ𝑡𝑢₁+Δ𝑡²2((𝐴+𝐵+𝐶)(𝑢₀−Δ𝑡3𝑢₁+Δ𝑡²12(𝐴+𝐵+𝐶)𝑢₀+))Δ𝑡²2𝑓(𝑡^𝑛)−Δ𝑡³6𝜕𝑓(𝑡^𝑛)+𝜕𝑡Δ𝑡⁴𝜕24²𝑓(𝑡^𝑛)𝜕𝑡²+𝑂(Δ𝑡⁵).(4.17)

Remark 4.1. By solving the two or three splitting steps, it is important to mention that each solution ^0𝑥002𝑑𝑢^{𝑐0𝑥002𝑑𝑐}, ^0𝑥0𝑢^02𝑑𝑐, and 𝑢 is corrected only once by using the boundary conditions.
Otherwise, an “overdoing” of the boundary conditions takes place.

4.3. 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 𝑢^𝑛+1,0−2𝑢^𝑛+𝑢^𝑛−1=Δ𝑡²(𝐴+𝐵)𝑢^𝑛,𝑢^𝑛+1,1−𝑢^𝑛+1,0=Δ𝑡²𝜂𝐴(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1𝑢),_𝑛+1−𝑢^𝑛+1,1=Δ𝑡²𝜂𝐵(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1),(4.18) where 𝜂∈(0.0,0.5) and 𝐴,𝐵 are the spatial discretized operators.

If we eliminate the intermediate values in (4.18), we obtain 𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1=Δ𝑡²(𝐴+𝐵)(𝜂𝑢^𝑛+1+(1−2𝜂)𝑢^𝑛+𝜂𝑢^𝑛−1)+𝐵_𝜂(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1),(4.19) where 𝐵_𝜂=𝜂²Δ𝑡⁴(𝐴𝐵) and thus 𝐵_𝜂(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1)=𝑂(Δ𝑡⁶).

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.

4.4. 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 𝑂(ℎ^𝑝), 𝑝=2,4, for the discretization in space, where ℎ=ℎ_𝑥=ℎ_𝑦 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 𝑢_𝑡𝑡−𝐴𝑢−(𝜕_𝑡𝑡𝑢−𝐴^{0𝑥002𝑑𝑐}𝑢)=𝑂(Δ𝑡²),(4.20)where 𝜕_𝑡𝑡 is a second-order discretization in time and ^0𝑥0𝐴^02𝑑𝑐 is the discretized fourth-order spatial operator.

Proof. We add (4.18) and obtain the following, see also [1]:𝜕_𝑡𝑡𝑢^𝑛−𝐴^{0𝑥002𝑑𝑐}(𝜃𝑢^𝑛+1+(1−2𝜃)𝑢^𝑛+𝜃𝑢^𝑛−1)−𝐵^{0𝑥002𝑑𝑐}(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1)=0,(4.21) where ^0𝑥0𝐵^02𝑑𝑐=𝜃²Δ𝑡𝐴^{20𝑥002𝑑𝑐}𝐴₁^{0𝑥002𝑑𝑐}₂.
Therefore, we obtain a splitting error of ^0𝑥0𝐵^02𝑑𝑐(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1).
Sufficient smoothness assumed that we have (𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1)=𝑂(Δ𝑡²), and we obtain ^0𝑥0𝐵^02𝑑𝑐(𝑢^𝑛+1−2𝑢^𝑛+𝑢^𝑛−1)=𝑂(Δ𝑡⁴).
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.

Theorem 4.4. The stability of the method is given by ‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁺_𝑡𝑢^𝑛‖²+𝒫⁺(𝑢^𝑛,𝜃)≤‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁺_𝑡𝑢⁰‖²+𝒫⁺(𝑢⁰,𝜃),(4.22) where 𝜃∈[0.25,0.5] and 𝒫^±(𝑢^𝑗,𝜃)∶=𝜃(𝐴^{0𝑥002𝑑𝑐}𝑢^𝑗,𝑢^𝑗)+𝜃(𝐴^{0𝑥002𝑑𝑐}𝑢^𝑗±1,𝑢^𝑗±1)+(1−2𝜃)(𝐴^{0𝑥002𝑑𝑐}𝑢^𝑗,𝑢^𝑗±1).

Proof. We have to proof the theorem for a test function 𝜕_𝑡𝑢^𝑛, where 𝜕_𝑡 denotes the central difference.
For 𝑗≥1, we have((1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})𝜕_𝑡𝑡𝑢^𝑗,𝜕_𝑡𝑢^𝑗)+(𝐴^{0𝑥002𝑑𝑐}(𝜃𝑢^𝑗+1−(1−2𝜃)𝑢^𝑗+𝜃𝑢^𝑗−1),𝜕_𝑡𝑢^𝑗)=0.(4.23)
Multiplying with Δ𝑡 and summarizing over 𝑗 yield^𝑛_𝑗=1((1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})𝜕_𝑡𝑡𝑢^𝑗,𝜕_𝑡𝑢^𝑗)Δ𝑡+(𝐴^{0𝑥002𝑑𝑐}((𝜃𝑢^𝑗+1−(1−2𝜃)𝑢^𝑗+𝜃𝑢^𝑗−1),𝜕_𝑡𝑢^𝑗),𝜕_𝑡𝑢^𝑗)Δ𝑡=0.(4.24) We can derive the identities ((1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})𝜕_𝑡𝑡𝑢^𝑗,𝜕_𝑡𝑢^𝑗1)Δ𝑡=2‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁺_𝑡𝑢^𝑗‖²−12‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁻_𝑡𝑢^ℎ‖²,(𝐴^{0𝑥002𝑑𝑐}(𝜃𝑢^𝑗+1−(1−2𝜃)𝑢^𝑗+𝜃𝑢^𝑗−1),𝜕_𝑡𝑢^𝑗1)Δ𝑡=2(𝒫⁺(𝑢^𝑗,𝜃)−𝒫⁻(𝑢^𝑗,𝜃)),(4.25) and obtain the result ‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁺_𝑡𝑢^𝑛‖²+𝒫⁺(𝑢^𝑛,𝜃)≤‖(1−Δ𝑡𝐵^{20𝑥002𝑑𝑐})^1/2𝜕⁺_𝑡𝑢⁰‖²+𝒫⁺(𝑢⁰,𝜃),(4.26) see also the idea of [1].

Remark 4.5. For 𝜃=1/12, we obtain a fourth-order method.
To compute the error of the local splitting, we have to use the multiplier 𝐴𝐴^{0𝑥002𝑑𝑐}₁^{0𝑥002𝑑𝑐}₂, 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].
(2) The truncation error is 𝑂(Δ𝑡²+ℎ^𝑝), 𝑝≥2, for 𝜃∈[0,0.5].
(3) 𝜃=1/12, we have a fourth-order method in time 𝑂(Δ𝑡²+ℎ^𝑝),𝑝≥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 𝐶𝐹𝐿=Δ𝑡²/Δ𝑥²_max𝐷_max, where 𝑥_max are the maximal spatial step and 𝐷_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.

4.5. 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 ℳ₄ into three parts: ℳ_𝑥𝑥, ℳ_𝑦𝑦, and ℳ_𝑥𝑦, where we have ℳ_𝑥𝑥=ℎ(𝜆+2𝜇)(1−²𝐷12^𝑥2)𝐷^𝑥20ℎ0𝜇(1−²𝐷12^𝑥2)𝐷^𝑥2,ℳ_𝑦𝑦=ℎ𝜇(1−²𝐷12^𝑦2)𝐷^𝑦20ℎ0(𝜆+2𝜇)(1−²𝐷12^𝑦2)𝐷^𝑦2,ℳ_𝑥𝑦=ℳ₄−ℳ_𝑥𝑥−ℳ_𝑦𝑦.(4.27) Our proposed split method has the following steps: 𝐔(1)𝜌^∗𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘Δ𝑡²=ℳ₄𝐔^𝑛_𝑗,𝑘+𝜃𝐟^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐟^𝑛_𝑗,𝑘+𝜃𝐟^𝑛−1_𝑗,𝑘,𝐔(2)𝜌^{∗∗𝑗,𝑘}−𝐔^∗_𝑗,𝑘Δ𝑡²=𝜃ℳ_𝑥𝑥(𝐔^∗∗_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘𝜃)+2ℳ_𝑥𝑦(𝐔^∗_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘𝐔),(3)𝜌^{𝑛+1𝑗,𝑘}−𝐔^∗∗_𝑗,𝑘Δ𝑡²=𝜃ℳ_𝑦𝑦(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘𝜃)+2ℳ_𝑥𝑦(𝐔^∗∗_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘).(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.

4.6. 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 𝐔(𝐱,𝑡) 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 𝐔 is 𝒪(Δ𝑡⁴), where it is assumed that Δ𝑡=𝒪(ℎ).

Proof. We assume in the following that 𝐟=(0,0)^𝑇. We add (4.28) and obtain, like in the scalar case [1], the following result for the discretized equations 𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘)−𝒩_4,𝜃(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)=𝟎,(4.30) where 𝒩_4,𝜃=𝜃²Δ𝑡²(ℳ_𝑥𝑥ℳ_𝑦𝑦+ℳ_𝑥𝑥ℳ_𝑥𝑦+ℳ_𝑥𝑦ℳ_𝑦𝑦)+𝜃³Δ𝑡⁴ℳ_𝑥𝑥ℳ_𝑦𝑦ℳ_𝑥𝑦.
We, therefore, obtain a splitting error of 𝒩_4,𝜃(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘).
For sufficient smoothness, we have (𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)=𝒪(Δ𝑡²) and we obtain 𝒩_4,𝜃(𝐔^𝑛+1−2𝐔^𝑛+𝐔^𝑛−1)=𝒪(Δ𝑡⁴).
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 𝐿₂(Ω).

In the following, we introduce the notation of the norms.

Remark 4.8. For our system, we extend the 𝐿₂-norm as ‖𝐔‖²_𝐿2=(𝐔,𝐔)_𝐿2=_Ω𝐔²𝑑𝑥,(4.31) where 𝐔²=𝑢²+𝑣² or 𝐔²=𝑢²+𝑣²+𝑤² in two and three dimensions.

Remark 4.9. For a discrete grid function 𝐔^𝑛_𝑗,𝑘, the 𝐿₂-norm is given as_Ω(𝐔^𝑛_𝑗𝑘)²𝑑𝑥=Δ𝑥^2𝑀_𝑗,𝑘𝑈^𝑛_𝑗,𝑘,(4.32) where Δ𝑥 is the uniform grid length in 𝑥 and 𝑦, 𝑀 is the number of grid points in 𝑥 and 𝑦. Further, 𝑈^𝑛_𝑗,𝑘 is the solution at grid point (𝑗,𝑘) and at time 𝑡^𝑛.

Remark 4.10. The matrix 𝒩_4,𝜃=𝜃²Δ𝑡²(ℳ_𝑥𝑥ℳ_𝑦𝑦+ℳ_𝑥𝑥ℳ_𝑥𝑦+ℳ_𝑥𝑦ℳ_𝑦𝑦)+𝜃³Δ𝑡⁴ℳ_𝑥𝑥ℳ_𝑦𝑦ℳ_𝑥𝑦,(4.33) where ℳ_𝑥𝑥, ℳ_𝑦𝑦, and ℳ_𝑥𝑦 are symmetrical and positive-definite matrices, therefore, the matrix 𝒩_4,𝜃 is also symmetrical and positive-definite.

Furthermore, we can estimate the norms and define a weighted norm, see [9, 10].

Remark 4.11. The energy norm is given as (𝒩_4,𝜃𝐔,𝐔)_𝐿2=_Ω(𝒩_4,𝜃𝐔𝐔)𝑑𝑥.(4.34) Consequently, we can denote ‖𝒩_4,𝜃𝐔‖≤𝜔‖𝐔‖∀𝐔∈𝐻^𝑑,(4.35) where 𝜔∈ℝ⁺ is the weight and 𝒩_4,𝜃 is bounded. 𝑑 is the dimension, and 𝐻 is Sobolev-space, see [11].

The stability of the fourth-order splitting method is given in the following theorem.

Theorem 4.12. Let 𝜃∈[0.25,0.5], then the implicit time-stepping algorithm, see (3.5), and the split procedure, see (4.28), are unconditionally stable. One can estimate the split procedure iteratively as ‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₊𝐔^𝑛_𝑗,𝑘‖²+𝒫⁺(𝐔^𝑛_𝑗,𝑘,𝜃)≤‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₊𝐔⁰_𝑗,𝑘‖²+𝒫⁺(𝐔⁰_𝑗,𝑘,𝜃),(4.36) where one has 𝒫^±(𝐔^𝑛_𝑗,𝑘,𝜃)∶=𝜃(ℳ₄𝐔^𝑛_𝑗,𝑘,𝐔^𝑛_𝑗,𝑘)+𝜃(ℳ₄𝐔^𝑛±1_𝑗,𝑘,𝐔^𝑛±1_𝑗,𝑘)+(1−2𝜃)(ℳ₄𝐔^𝑛_𝑗,𝑘,𝐔^𝑛±1_𝑗,𝑘) and 𝒫^±≥0 for 𝜃∈[0.25,0.5]. Further, 1−Δ𝑡𝜔^{20𝑥002𝑑𝑐}∈ℝ⁺ is the factor for the weighted norm (ℐ−Δ𝑡²𝒩_4,𝜃)𝐔≤𝜔^{0𝑥002𝑑𝑐}𝐔 for all 𝐔∈𝐻^𝑑.

We have to prove the iterative estimate for the split procedure and the proof is given as follows.

Proof. To obtain an energy estimate for the scheme, we multiply with a test-function 𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘.
The following result is given for the discretized equations, see also (4.30):(ℐ−Δ𝑡²𝒩_4,𝜃)𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘)=𝟎.(4.37)
So for 𝑛≥1, we can rewrite (4.37) for the stability proof as follows: ((ℐ−Δ𝑡²𝒩_4,𝜃)𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘,𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘)−(ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘),𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘)=𝟎.(4.38)
Multiplying with Δ𝑡 and summarizing over the time levels, we obtain _𝑛((ℐ−Δ𝑡²𝒩_4,𝜃)𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘,𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘)Δ𝑡−_𝑛(ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘),𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘)Δ𝑡=𝟎,(4.39) for each term of the sum, one can derive the following identities. So for ℐ−Δ𝑡²𝒩_4,𝜃, we have ((ℐ−Δ𝑡²𝒩_4,𝜃)𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘,𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘1)Δ𝑡=2((ℐ−Δ𝑡²𝒩_4,𝜃)(𝐷^𝑡₊−𝐷^𝑡₋)𝐔^𝑛_𝑗,𝑘,(𝐷^𝑡₊+𝐷^𝑡₋)𝐔^𝑛_𝑗,𝑘)=_Ω((ℐ−Δ𝑡²𝒩_4,𝜃)(𝐷^𝑡₊−𝐷^𝑡₋))^𝑇(𝐷^𝑡₊+𝐷^𝑡₋)𝐔^𝑛_𝑗,𝑘𝑑𝑥≤(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})_Ω(𝐷^𝑡₊𝐔^𝑛_𝑗,𝑘)²(𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘)²𝑑𝑥,(4.40) where the operator ℐ−Δ𝑡²𝒩_4,𝜃 is symmetric and positive-definite and we can apply the weighted norm, see Remark 4.11 and [11].
We obtain the following result: (1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})_Ω(𝐷^𝑡₊𝐔^𝑛_𝑗,𝑘)²(𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘)²1𝑑𝑥=2‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₊𝐔^𝑛_𝑗,𝑘‖²−12‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘‖².(4.41)
Further, for −ℳ₄, we have (−ℳ₄(𝜃𝐔^𝑛+1_𝑗,𝑘+(1−2𝜃)𝐔^𝑛_𝑗,𝑘+𝜃𝐔^𝑛−1_𝑗,𝑘),𝐷^𝑡₀𝐔^𝑛_𝑗,𝑘1)Δ𝑡=2(𝒫⁺(𝐔^𝑛_𝑗,𝑘,𝜃)−𝒫⁻(𝐔^𝑛_𝑗,𝑘,𝜃)).(4.42)
Due to the result of the operators: 𝒫⁻(𝐔^𝑛_𝑗,𝑘,𝜃)=𝒫⁺(𝐔^𝑛−1_𝑗,𝑘,𝜃),𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘=𝐷^𝑡₊𝐔^𝑛−1_𝑗,𝑘,(4.43) we can recursively derive the following result: ‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₊𝐔^𝑛_𝑗,𝑘‖²+𝒫⁺(𝐔^𝑛_𝑗,𝑘,𝜃)≤‖(1−Δ𝑡𝜔^{20𝑥002𝑑𝑐})^1/2𝐷^𝑡₊𝐔⁰_𝑗,𝑘‖²+𝒫⁺(𝐔⁰_𝑗,𝑘,𝜃),(4.44) where for 𝜃∈[0.25,0.5], we have 𝒫⁺(𝐔^𝑛_𝑗,𝑘,𝜃)≥0 for all 𝑛∈ℕ⁺, and, therefore, we have the unconditional stability. The scalar proof is also presented in the work of [1].

Remark 4.13. For 𝜃=1/12, the split method is fourth-order accurate in time and space.

See the following theorem.

Theorem 4.14. One obtains a fourth-order accurate scheme in time and space for the split method, see (4.28), when 𝜃=1/12. That reads 𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−1ℳ12⁴(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)+ℳ₄𝐔^𝑛_𝑗,𝑘+𝒩_4,𝜃(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)=𝟎,(4.45) where ℳ₄ is a fourth-order discretization scheme in space.

The proof is given as follows.

Proof. We consider the following Taylor-expansion: 𝜕_𝑡𝑡𝐔^𝑛_𝑗,𝑘=𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−Δ𝑡²𝜕12^{𝑡𝑡𝑡𝑡}𝐔^𝑛_𝑗,𝑘+𝒪(Δ𝑡⁴).(4.46) Furthermore, we have 𝜕_{𝑡𝑡𝑡𝑡}𝐔^𝑛_𝑗,𝑘≈ℳ₄𝜕_𝑡𝑡𝐔^𝑛_𝑗,𝑘,(4.47) and we can rewrite (4.46) as 𝜕_𝑡𝑡𝐔^𝑛_𝑗,𝑘≈𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−Δ𝑡²ℳ12⁴𝜕_𝑡𝑡𝐔^𝑛_𝑗,𝑘+𝒪(Δ𝑡⁴)≈𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−Δ𝑡²ℳ12⁴(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)+𝒪(Δ𝑡⁴).(4.48)
So the fourth-order time-stepping algorithm can be formulated as 𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−1ℳ12⁴(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)−ℳ₄𝐔^𝑛_𝑗,𝑘=𝟎.(4.49)
The split method (4.28) becomes 𝐷^𝑡₊𝐷^𝑡₋𝐔^𝑛_𝑗,𝑘−1ℳ12⁴(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)−ℳ₄𝐔^𝑛_𝑗,𝑘−𝒩_4,1/12(𝐔^𝑛+1_𝑗,𝑘−2𝐔^𝑛_𝑗,𝑘+𝐔^𝑛−1_𝑗,𝑘)=𝟎,(4.50) and we obtain a fourth-order split scheme (cf. the scalar case [1]).

Remark 4.15. As follows form Theorem 4.14, we obtain a fourth-order in time for 𝜃=1/12. For the stability analysis, the method is conditional stable for 𝜃∈(0,0.25). So the splitting method will not restrict our stability condition for the fourth-order method with 𝜃=1/12.

Our theoretical results are verified by the following numerical examples.

5. Numerical Experiments

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

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

5.1.1. 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: 𝜕_𝑡𝑡𝑢=𝐷²₁𝜕_𝑥𝑥𝑢+𝐷²₂𝜕_𝑦𝑦𝑢,inΩ×(0,𝑇),𝑢(𝑥,𝑦,0)=𝑢₀1(𝑥,𝑦)=sin(𝐷₁1𝜋𝑥)sin(𝐷₂𝜕𝜋𝑦),onΩ,_𝑡𝑢(𝑥,𝑦,0)=𝑢₁1(𝑥,𝑦)=0,onΩ,𝑢(𝑥,𝑦,𝑡)=sin(𝐷₁1𝜋𝑥)sin(𝐷₂𝜋𝑦)cos(√2𝜋𝑡),on𝜕Ω×(0,𝑇),(5.1) where the initial conditions can be written as 𝑢(𝑥,𝑦,𝑡^𝑛)=𝑢₀(𝑥,𝑦) and 𝑢(𝑥,𝑦,𝑡^𝑛−1)=𝑢(𝑥,𝑦,𝑡^𝑛+1)=𝑢(𝑥,𝑦,Δ𝑡).

The analytical solution is given by 𝑢_ana1(𝑥,𝑦,𝑡)=sin(𝐷₁1𝜋𝑥)sin(𝐷₂𝜋𝑦)cos(√2𝜋𝑡).(5.2)

For the approximation error, we choose the 𝐿₁-norm.

The 𝐿₁-norm is given by err_𝐿1∶=_{𝑖,𝑗=1,…,𝑚}𝑉_𝑖,𝑗|𝑢(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛)−𝑢_ana(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛)|,(5.3) where 𝑢(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛) is the numerical and 𝑢_ana(𝑥_𝑖,𝑦_𝑗,𝑡^𝑛) is the analytical solution and 𝑉_𝑖,𝑗=Δ𝑥Δ𝑦.

Our test examples are organized as follows. (1)The non-stiff case. We choose 𝐷₁=𝐷₂=1 with a rectangle as our model domain Ω=[0,1]×[0,1]. We discretize with Δ𝑥=1/16 and Δ𝑦=1/16 and Δ𝑡=1/32 and choose our parameter 𝜂 between 0≤𝜂≤1. The exemplary function values 𝑢_num and 𝑢_ana are taken from the center of our domain.(2)The stiff case. We choose 𝐷₁=𝐷₂=0.01 with a rectangle as our model domain Ω=[0,1]×[0,1]. We discretize with Δ𝑥=1/32 and Δ𝑦=1/32 and Δ𝑡=1/64 and choose our parameter 𝜂 between 0≤𝜂≤1. The exemplary function values 𝑢_num and 𝑢_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.

Figure 1 

Numerical error for standard and modified methods, with respect to the 𝜂 parameter and Dirichlet boundary conditions.

(a)

(b)

Figure 2 

Numerical resolution of the wave equation: numerical approximation (a) and error functions (b) for the Dirichlet boundary conditions (Δ𝑥=Δ𝑦=1/32, Δ𝑡=1/64, 𝐷1=1,𝐷2=1) (classical method).

The numerical errors for the stiff case with Dirichlet boundary conditions are presented in Figure 3 and their results in Figure 4.

Figure 3 

Numerical error for standard and modified methods, with respect to the 𝜂 parameter, Dirichlet boundary conditions and stiff case (Δ𝑥=Δ𝑦=1/32, Δ𝑡=1/64, 𝐷1=1,𝐷2=0.01).

(a)

(b)

Figure 4 

Numerical approximation and error function for the Dirichlet boundary in the stiff case (Δ𝑥=Δ𝑦=1/32, Δ𝑡=1/64, 𝐷1=1,𝐷2=0.01.).

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.

In the next test example, we study the Neumann boundary conditions.

5.1.2. 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: 𝜕_𝑡𝑡𝑢=𝐷²₁𝜕_𝑥𝑥𝑢+𝐷²₂𝜕_𝑦𝑦𝑢,inΩ×(0,𝑇),𝑢(𝑥,𝑦,0)=𝑢₀1(𝑥,𝑦)=sin(𝐷₁1𝜋𝑥)sin(𝐷₂𝜕𝜋𝑦),onΩ,_𝑡𝑢(𝑥,𝑦,0)=𝑢₁(𝑥,𝑦)=0,onΩ,𝜕𝑢(𝑥,𝑦,𝑡)=𝜕𝑛𝜕𝑢_analy(𝑥,𝑦,𝑡)𝜕𝑛=0,on𝜕Ω×(0,𝑇),(5.4) where Ω=[0,1]×[0,1]. 𝐷₁=1, 𝐷₂=0.5 and the initial conditions can be written as 𝑢(𝑥,𝑦,𝑡^𝑛)=𝑢₀(𝑥,𝑦) and 𝑢(𝑥,𝑦,𝑡^𝑛−1)=𝑢(𝑥,𝑦,𝑡^𝑛+1)=𝑢(𝑥,𝑦,Δ𝑡).

The analytical solution is given as𝑐_analy1(𝑥,𝑦,𝑡)=sin(𝐷₁1𝜋𝑥)sin(𝐷₂𝜋𝑦)cos(√2𝜋𝑡).(5.5)

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.

Figure 5 

Numerical error for standard and modified methods, with respect to the 𝜂 parameter and Neumann boundary conditions.

(a)

(b)

Figure 6 

Numerical resolution of the wave equation: numerical approximation (a) and error functions (b) for the Neumann boundary condition (Δ𝑥=Δ𝑦=1/32, Δ𝑡=1/64, 𝐷1=1,𝐷2=1) (classical method).

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.

5.1.3. Test Example 3: Spatial-Dependent Wave Equation

In this experiment, we apply our method to the spatial-dependent problem, given by 𝜕_𝑡𝑡𝑢=𝐷₁(𝑥,𝑦)𝜕_𝑥𝑥𝑢+𝐷₂(𝑥,𝑦)𝜕_𝑦𝑦𝑢,inΩ×(0,𝑇),𝑢(𝑥,𝑦,𝑡^𝑛)=𝑢₀,𝜕_𝑡𝑢(𝑥,𝑦,𝑡^𝑛)=𝑢₁,on𝜕Ω×(0,𝑇),𝑢(𝑥,𝑦,𝑡)=𝑢₂,on𝜕Ω×(0,𝑇),(5.6) where 𝐷₁(𝑥,𝑦)=0.1𝑥+0.01𝑦+0.01, 𝐷₂(𝑥,𝑦)=0.01𝑥+0.1𝑦+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.

Remark 5.3. We have assumed the following CFL condition: Δ𝑡<0.5min(Δ𝑥,Δ𝑦)max_{𝑥,𝑦∈Ω}(𝐷₁(𝑥,𝑦),𝐷₂(𝑥,𝑦)).(5.7)

For the test example, we define our model domain as a rectangle Ω=[0,1]×[0,1].

The reference solution is obtained by executing the finite differences method and setting Δ𝑥=1/256,Δ𝑦=1/256, and the time step Δ𝑡=1/256<0.390625.

The model domain is given by a rectangle with Δ𝑥=1/16 and Δ𝑦=1/32. The time steps are given by Δ𝑡=1/16 and 0≤𝜂≤1.

The numerical errors for the spatial-dependent parameters with Dirichlet boundary conditions are presented in Figure 7 and their results in Figure 8.

Figure 7 

Numerical error for standard and modified methods, with respect to the 𝜂 parameter, spatial-dependent parameters, and Dirichlet boundary conditions.

(a)

(b)

Figure 8 

Dirichlet boundary condition: numerical solution and error function for the spatial-dependent test example.

The numerical errors for the spatial-dependent parameters with Neumann boundary conditions are presented in Figure 9 and their results in Figure 10.

Figure 9 

Numerical error for standard and modified methods with respect to the 𝜂 parameter, spatial-dependent parameters, and Neumann boundary conditions.

(a)

(b)

Figure 10 

Neumann boundary condition: numerical solution and error function for the spatial-dependent test example.

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.

5.2. 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 𝐟 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.

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.

5.2.1. 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 𝑡=0, we know the value at level 𝑛=0 as 𝐔⁰=𝐠₀. The value at level 𝑛=−1 can be obtained by Taylor expansion as 𝐔⁻¹=𝐔⁰−Δ𝑡𝜕_𝑡𝐔⁰+Δ𝑡²2𝜕_𝑡𝑡𝐔⁰−Δ𝑡³6𝜕_𝑡𝑡𝑡𝐔⁰+Δ𝑡⁴𝜕24^{𝑡𝑡𝑡𝑡}𝐔⁰+𝑂(Δ𝑡⁵),(5.8) where we use 𝜕_𝑡𝐔⁰_𝑗,𝑘=𝐠_1𝑗,𝑘𝜕,(5.9a)_𝑡𝑡𝐔⁰_𝑗,𝑘≈1𝜌(ℳ₄𝐠_0𝑗,𝑘)+𝐟_𝑗,𝑘𝜕),(5.9b)_𝑡𝑡𝑡𝐔⁰_𝑗,𝑘≈1𝜌(ℳ₄𝐠_1𝑗,𝑘)+𝜕_𝑡𝐟⁰_𝑗,𝑘𝜕),(5.9c)_{𝑡𝑡𝑡𝑡}𝐔⁰_𝑗,𝑘≈1𝜌(ℳ²₂𝐠_0𝑗,𝑘)+ℳ₄𝐟⁰_𝑗,𝑘+𝜕_𝑡𝑡𝐟⁰_𝑗,𝑘),(5.9d) and also for (5.9c) and (5.9d), 𝜕_𝑡𝐟⁰_𝑗,𝑘≈𝐟¹_𝑗,𝑘−𝐟⁻¹_𝑗,𝑘𝜕2Δ𝑡,(5.9e)_𝑡𝑡𝐟⁰_𝑗,𝑘≈𝐟¹_𝑗,𝑘−2𝐟⁰_𝑗,𝑘+𝐟⁻¹_𝑗,𝑘Δ𝑡².(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 𝐔^∗, 𝐔^∗∗, and for 𝐔^𝑛+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.2.2. Test Example

For the first test case, we use a forcing function 𝐟=(sin(𝑡−𝑥)sin(𝑦)−2𝜇sin(𝑡−𝑥)sin(𝑦)−(𝜆+𝜇)(cos(𝑥)cos(𝑡−𝑦)+sin(𝑡−𝑥)sin(𝑦)),sin(𝑡−𝑦)sin(𝑥)−2𝑉𝑠²sin(𝑥)sin(𝑡−𝑦)−(𝜆+𝜇)(cos(𝑡−𝑥)cos(𝑦)+sin(𝑦)sin(𝑡−𝑦)))^𝑇,(5.10) giving the analytical solution 𝐔^true=(sin(𝑥−𝑡)sin(𝑦),sin(𝑦−𝑡)sin(𝑥))^𝑇.(5.11) Using the split method we solved (2.4a) on a domain 𝑥×𝑦=[−11]×[−11] up to 𝑡=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 ‖𝐔_𝑗,𝑘‖=max(max_𝑗,𝑘|𝑢_𝑗,𝑘|,max_𝑗,𝑘|𝑣_𝑗,𝑘|). As can be seen we get the expected 4th order convergence for problems with smooth solutions.

Table 1 

Errors in max-norm for decreasing ℎ and smooth analytical solution 𝐔^true. Convergence rate indicates a fourth-order convergence for the split scheme.

To check the influence of the splitting error 𝒩_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.

Table 2 

Errors in max-norm for decreasing ℎ and smooth analytical solution 𝐔^true and using the non-split scheme. Comparing with Table 1, we see that the splitting error is very small for this case.

5.2.3. Singular Forcing Terms

In seismology and acoustics it is common to use spatially singular forcing terms which can look like𝐟=𝐅𝛿(𝐱)𝑔(𝑡),(5.12)where 𝐅 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 (𝛿_ℎ²) and one 4th (𝛿_ℎ⁴) order regularized approximations of the one dimensional Dirac function, 𝛿_ℎ²(𝑥^{0𝑥002𝑑𝑐}1)=ℎ1+𝑥^{0𝑥002𝑑𝑐},−ℎ≤𝑥^{0𝑥002𝑑𝑐}<0,1−𝑥^{0𝑥002𝑑𝑐},0≤𝑥^{0𝑥002𝑑𝑐}𝛿<ℎ,0,elsewhere,(5.13)_ℎ4(𝑥^{0𝑥002𝑑𝑐}1)=ℎ1+116^{0𝑥002𝑑𝑐}𝑥+58^{0𝑥002𝑑𝑐}𝑥²+16^{0𝑥002𝑑𝑐}𝑥³,−2ℎ≤𝑥^{0𝑥002𝑑𝑐}1<−ℎ,1+2^{0𝑥002𝑑𝑐}𝑥−𝑥^{0𝑥002𝑑𝑐2}−12^{0𝑥002𝑑𝑐}𝑥³,−ℎ≤𝑥^{0𝑥002𝑑𝑐}1<0,1−2^{0𝑥002𝑑𝑐}𝑥−𝑥^{0𝑥002𝑑𝑐2}+12^{0𝑥002𝑑𝑐}𝑥³,0≤𝑥^{0𝑥002𝑑𝑐}<ℎ,1−116^{0𝑥002𝑑𝑐}𝑥+𝑥^{0𝑥002𝑑𝑐2}−16^{0𝑥002𝑑𝑐}𝑥³,ℎ≤𝑥^{0𝑥002𝑑𝑐}<2ℎ,0,elsewhere,(5.14)where in the above ^0𝑥0𝑥^02𝑑𝑐=𝑥/ℎ. The two and three dimensional Dirac functions are then approximated as 𝛿_ℎ^2,4(𝑥^{0𝑥002𝑑𝑐})𝛿_ℎ^2,4(𝑦^{0𝑥002𝑑𝑐}) and 𝛿_ℎ^2,4(𝑥^{0𝑥002𝑑𝑐})𝛿_ℎ^2,4(𝑦^{0𝑥002𝑑𝑐})𝛿_ℎ^2,4(𝑧^{0𝑥002𝑑𝑐}). The chosen time dependence was a smooth function given by 𝑔(𝑡)=exp(−1𝑡(1−𝑡)),0≤𝑡<1,0,elsewhere,(5.15) which is 𝐶^∞. 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 𝑡=0.8 for different grid sizes ℎ and the two approximations 𝛿_ℎ² and 𝛿_ℎ⁴. As can be seen the error is smooth and converges a small distance away from the source. However, using 𝛿_ℎ² limits the convergence to 2nd order, while using 𝛿_ℎ⁴ gives the full 4th order convergence away from the singular source. When 𝑡>1 the forcing goes to zero and the solution will be smooth everywhere. Table 3 shows the convergence behavior at time 𝑡=1.1 for four different grids. Note that the full convergence is achieved even if the lower order 𝛿_ℎ² is used as an approximation for the Dirac function. The convergence rate approaches 4 as we refine the grids, even though the solution was singular up to time 𝑡=1.

Table 3 

Errors in max-norm for decreasing ℎ and analytical solution 𝐔^true. Convergence rate approaches 4th order after the singular forcing term goes to zero.

Figure 11 

The 2-logarithm of the error along a line going through the source point for a point force located at 𝑥=0, 𝑦=0, and approximated in space by (5.14). Note that the error decays as 𝒪(ℎ4) away from the source, but not near it. The grid sizes were ℎ=0.05(−⋅),0.025(⋅),0.0125(−),0.006125(∗). The numerical quadrature had an absolute error of approximately 10−11≈2−36, so the error cannot be resolved beneath that limit.

Figure 12 

The 2-logarithm of the error along a line going through the source point for a point force located at 𝑥=0, 𝑦=0, and approximated in space by (5.13). Note that the error only decays as 𝒪(ℎ2) away from the source. The grid sizes were ℎ=0.05(−⋅),0.025(⋅),0.0125(−),0.006125(∗).

Remark 5.6. For a two dimensional problem the 4th order explicit method (3.10) can be implemented using approximately 160 floating point operations (flops) per grid point. For example, the split method requires approximately 120 flops (first step) plus 2 times 68 flops (second and third step) for a total of 256 flops. This increase of about 60% in the number of flops is somewhat offset by the larger time steps allowed by the split method, especially for smooth material properties, making the two methods roughly comparable in computational cost.

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

5.3.1. The Splitting Scheme

In three dimensions a 4th order difference approximation of the right hand side operator becomesℳ₄=ℎ(𝜆+2𝜇)(1−²𝐷12^𝑥2)𝐷^𝑥2ℎ+𝜇(1−²𝐷12^𝑦2𝐷^𝑦2ℎ+1−²𝐷12^𝑧2)𝐷^𝑧2ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑧2)𝐷^𝑧₀ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(𝜆+2𝜇)(1−²𝐷12^𝑦2)𝐷^𝑦2ℎ+𝜇(1−²𝐷12^𝑥2)𝐷^𝑥2ℎ+𝜇(1−²𝐷12^𝑧2)𝐷^𝑧2ℎ(𝜆+𝜇)(1−²6𝐷^𝑧2)𝐷^𝑧₀ℎ(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(𝜆+𝜇)(1−²6𝐷^𝑥2)𝐷^𝑥₀ℎ(1−²6𝐷^𝑧2)𝐷^𝑧₀ℎ(𝜆+𝜇)(1−²6𝐷^𝑦2)𝐷^𝑦₀ℎ(1−²6𝐷^𝑧2)𝐷^𝑧₀ℎ(𝜆+2𝜇)(1−²𝐷12^𝑧2𝐷^𝑧2ℎ+𝜇(1−²𝐷12^𝑥2)𝐷^𝑥2ℎ+1−²𝐷12^𝑦2)𝐷^𝑦2,(5.16) operating on grid functions 𝐔^𝑛_{𝑗,𝑘,𝑙} defined at grid points 𝑥_𝑗,𝑦_𝑘,𝑧_𝑙,𝑡_𝑛 similarly to the two dimensional case. We can split ℳ₄ into six parts; ℳ_𝑥𝑥, ℳ_𝑦𝑦, ℳ_𝑧𝑧 containing the three second order directional difference operators, and ℳ_𝑥𝑦, ℳ_𝑦𝑧, ℳ_𝑥𝑧 containing the mixed difference operators.

There are a number of different ways we could split this scheme, depending on how we treat the mixed derivative terms. We have chosen to implement the following split scheme in three dimensions: 𝐔(1)𝜌^{∗𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙}Δ𝑡²=ℳ₄𝐔^𝑛_{𝑗,𝑘,𝑙}+𝜃𝐟^𝑛+1_{𝑗,𝑘,𝑙}+(1−2𝜃)𝐟^𝑛_{𝑗,𝑘,𝑙}+𝜃𝐟^𝑛−1_{𝑗,𝑘,𝑙}𝐔(2)𝜌^{∗∗𝑗,𝑘,𝑙}−𝐔^∗_{𝑗,𝑘,𝑙}Δ𝑡²=𝜃ℳ_𝑥𝑥(𝐔^∗∗_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙}𝜃)+2(ℳ_𝑥𝑦+ℳ_𝑥𝑧)(𝐔^∗_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙})𝐔(3)𝜌^{∗∗∗𝑗,𝑘,𝑙}−𝐔^∗∗_{𝑗,𝑘,𝑙}Δ𝑡²=𝜃ℳ_𝑥𝑥(𝐔^∗∗∗_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙}𝜃)+2(ℳ_𝑥𝑦+ℳ_𝑦𝑧)(𝐔^∗∗_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙})𝐔(4)𝜌^{𝑛+1𝑗,𝑘,𝑙}−𝐔^∗∗∗_{𝑗,𝑘,𝑙}Δ𝑡²=𝜃ℳ_𝑥𝑥(𝐔^𝑛+1_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙}𝜃)+2(ℳ_𝑥𝑧+ℳ_𝑦𝑧)(𝐔^∗∗∗_{𝑗,𝑘,𝑙}−2𝐔^𝑛_{𝑗,𝑘,𝑙}+𝐔^𝑛−1_{𝑗,𝑘,𝑙}).(5.17)

The properties such as splitting error, accuracy, stability, and so forth, for the three dimensional case are similar to the two dimensional case treated in the earlier sections.

5.3.2. Testing the Three Dimensional Scheme

We have done some numerical experiments with the three dimensional scheme in order to test the convergence and stability. We used a forcing 𝐟=−(−1+𝜆+4𝜇)sin(𝑡−𝑥)sin(𝑦)sin(𝑧)−(𝜆+𝜇)cos(𝑥)(2sin(𝑡)sin(𝑦)sin(𝑧)+cos(𝑡)sin(𝑦+𝑧)),−(−1+𝜆+4𝜇)sin(𝑥)sin(𝑡−𝑦)sin(𝑧)−(𝜆+𝜇)cos(𝑦)(2sin(𝑡)sin(𝑥)sin(𝑧)+cos(𝑡)sin(𝑥+𝑧)),−(𝜆+𝜇)cos(𝑡−𝑦)cos(𝑧)sin(𝑥)−sin(𝑦)((𝜆+𝜇)cos(𝑡−𝑥)cos(𝑧)+(−1+𝜆+4𝜇)sin(𝑥)sin(𝑡−𝑧)))^𝑇,(5.18) giving the analytical solution 𝐔^true=(sin(𝑥−𝑡)sin(𝑦)sin(𝑧),sin(𝑦−𝑡)sin(𝑥)sin(𝑧),sin(𝑧−𝑡)sin(𝑥)sin(𝑦))^𝑇.(5.19) As earlier we tested this for a number of different grid sizes. Using the same two sets of material parameters as for the two dimensional case we ran up until 𝑡=2 and checked the max error for all components of the solution. The results are given in Table 4. We have also tested the three dimensional scheme using singular forcing functions approximated using (5.13) and (5.14). The results are very similar to the two dimensional case and we have therefore omitted them here.

Table 4 

Errors in max-norm for decreasing ℎ and smooth analytical solution 𝐔^true. Convergence rate indicates 4th order convergence for the three dimensional split scheme.

6. Conclusion

We have presented time splitting methods for the scalar and vectorial wave equation. The contributions of this article concerns the higher order splitting methods, based on LOD 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.