# Higher-Order Splitting Method

Authors: Jürgen Geiser
Journal: International Journal of Mathematics and Mathematical Sciences
 (2008)

---

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


---

## Body

## 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 Section2. 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
(2.1)∂ttu=D∇⋅∇u,inΩ,u(x,0)=u0(x),ut(x,0)=u1(x),inΩ.
The unknown function u=u(x,t) is considered to be in Ω×(0,T)⊂ℝd×ℝ, where the spatial dimension is given by d.For three dimensions, we define the diffusion tensor as(2.2)D=(D1000D2000D3)∈ℝ3,+×ℝ3,+,
which describes the wave
propagation. Further, the diffusion tensor D is given anisotropic, with D1,D2,D3∈ℝ+ for D1,D2≥D3.
The functions u0(x) and u1(x) are the initial conditions for the wave
equation.We deal with the following boundary conditions:(2.3)u(x,t)=u3,Dirichletboundarycondition,∂u(x,t)∂n=0,Neumannboundarycondition,D∇u(x,t)=uout,Outflowboundarycondition,
where all boundary conditions are on ∂Ω×T.


### 2.2. Elastic Wave Propagation

We consider the initial-value problem for the elastic
wave equation for constant coefficients, given  as(2.4a)ρ∂ttU=μ∇2U+(λ+μ)∇(∇⋅U)+f,(2.4b)U(t=0,x)=g0(x),(2.4c)∂tU(t=0,x)=g1(x),
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,g0,
and g1 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(2.5)f=Fδ(x)g(t),
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.


## 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
(2.1)∂ttu=D∇⋅∇u,inΩ,u(x,0)=u0(x),ut(x,0)=u1(x),inΩ.
The unknown function u=u(x,t) is considered to be in Ω×(0,T)⊂ℝd×ℝ, where the spatial dimension is given by d.For three dimensions, we define the diffusion tensor as(2.2)D=(D1000D2000D3)∈ℝ3,+×ℝ3,+,
which describes the wave
propagation. Further, the diffusion tensor D is given anisotropic, with D1,D2,D3∈ℝ+ for D1,D2≥D3.
The functions u0(x) and u1(x) are the initial conditions for the wave
equation.We deal with the following boundary conditions:(2.3)u(x,t)=u3,Dirichletboundarycondition,∂u(x,t)∂n=0,Neumannboundarycondition,D∇u(x,t)=uout,Outflowboundarycondition,
where all boundary conditions are on ∂Ω×T.


## 2.2. Elastic Wave Propagation

We consider the initial-value problem for the elastic
wave equation for constant coefficients, given  as(2.4a)ρ∂ttU=μ∇2U+(λ+μ)∇(∇⋅U)+f,(2.4b)U(t=0,x)=g0(x),(2.4c)∂tU(t=0,x)=g1(x),
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,g0,
and g1 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(2.5)f=Fδ(x)g(t),
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.


## 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(3.1)Utt,i=Uin+1−2Uin+Uin−1Δt2,U(tn)=u(x,t),Ut(tn)=ut(x,t),
where the index i refers to the space point xi and Δt=tn+1−tn is the time step. We apply finite difference
methods for the spatial discretization. The spatial terms and the initial
conditions are given as
(3.2)Uxx,n=Ui+1n−2Uin+Ui−1nΔx2,U(tn)=u(x,t),Ut(tn)=ut(x,t),
where the index n refers to the time tn and Δx=xi+1−xi is the grid width.Then the two-dimensional equation,(3.3)utt=D1uxx+D2uyyinΩ,u(x,y,0)=u0(x,y),ut(x,y,0)=u1(x,y),u(x,y,t)=u2on∂Ω,
is discretized with the unconditionally stable implicit η-method, see [5],
(3.4)Ui,jn+1−2Ui,jn+Ui,jn−1Δt2=D1Δx2(η(Ui+1,jn+1−2Ui,jn+1+Ui−1,jn+1)+(1−2η)(Ui+1,jn−2Ui,jn+Ui−1,jn)+η(Ui+1,jn−1−2Ui,jn−1+Ui−1,jn−1))+D2Δy2(η(Ui,j+1n−2Ui,jn+Ui,j−1n)+(1−2η)(Ui,j+1n−2Ui,jn+Ui,j−1n)+η(Ui,j+1n−1−2Ui,jn−1+Ui,j−1n−1)),where Δx and Δy are the grid width in x and y and 0≤η≤1.
The initial conditions are given by U(x,y,tn)=u(x,y,tn) and U(x,y,tn−1)=u(x,y,tn)−Δtut(x,y,tn).These discretization schemes are adopted to the
operator splitting schemes.On the finite differences grid,Δt corresponds to the time step, and hx,hy are the grid sizes in the different spatial
directions. The time nΔt is denoted by tn,
and i,j refer to the spatial coordinates of the grid
point (ihx,jhy).
Let un denote the grid function on the time level n,
and ui,jn be the specific value of un at point i,j.In Section3.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 pointsxj=jh,yk=kh,tn=nΔt, where h>0 is the spatial grid size and Δt the time step. Defining the grid function Uj,kn=U(xj,yk,tn),
the basic explicit scheme is(3.5)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ2Uj,kn+fj,kn,where ℳ2 is a difference operator(3.6)ℳ2=((λ+2μ)Dx2+μDy2(λ+μ)D0xD0y(λ+μ)D0xD0y(λ+2μ)Dy2+μDx2),
and we use the standard
difference operator notation:
(3.7)D+xvj,k=1h(vj+1,k−vj,k),D−xvj,k=D+xvj−1,k,D0x=12(D+x+D−x),Dx2=D+xD−xℳ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](3.8)Δt<hλ+3μ.Replacingℳ2 with ℳ4,
a fourth-order difference operator given by(3.9)ℳ4=((λ+2μ)(1−h212Dx2)Dx2+μ(1−h212Dy2)Dy2(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+2μ)(1−h212Dy2)Dy2+μ(1−h212Dx2)Dx2),
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
(3.10)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ4Uj,kn+fj,kn+Δt212(ℳ22Uj,kn+ℳ2fi,jn+∂ttfi,jn),where ℳ22 is a second-order approximation to the squared
right-hand side operator in (2.4a). As it only needs to be second-order accurate, ℳ22 has the same extent in space as ℳ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:(3.11)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)+θfj,kn+1+(1−2θ)fj,kn+θfj,kn−1.
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.


## 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(3.1)Utt,i=Uin+1−2Uin+Uin−1Δt2,U(tn)=u(x,t),Ut(tn)=ut(x,t),
where the index i refers to the space point xi and Δt=tn+1−tn is the time step. We apply finite difference
methods for the spatial discretization. The spatial terms and the initial
conditions are given as
(3.2)Uxx,n=Ui+1n−2Uin+Ui−1nΔx2,U(tn)=u(x,t),Ut(tn)=ut(x,t),
where the index n refers to the time tn and Δx=xi+1−xi is the grid width.Then the two-dimensional equation,(3.3)utt=D1uxx+D2uyyinΩ,u(x,y,0)=u0(x,y),ut(x,y,0)=u1(x,y),u(x,y,t)=u2on∂Ω,
is discretized with the unconditionally stable implicit η-method, see [5],
(3.4)Ui,jn+1−2Ui,jn+Ui,jn−1Δt2=D1Δx2(η(Ui+1,jn+1−2Ui,jn+1+Ui−1,jn+1)+(1−2η)(Ui+1,jn−2Ui,jn+Ui−1,jn)+η(Ui+1,jn−1−2Ui,jn−1+Ui−1,jn−1))+D2Δy2(η(Ui,j+1n−2Ui,jn+Ui,j−1n)+(1−2η)(Ui,j+1n−2Ui,jn+Ui,j−1n)+η(Ui,j+1n−1−2Ui,jn−1+Ui,j−1n−1)),where Δx and Δy are the grid width in x and y and 0≤η≤1.
The initial conditions are given by U(x,y,tn)=u(x,y,tn) and U(x,y,tn−1)=u(x,y,tn)−Δtut(x,y,tn).These discretization schemes are adopted to the
operator splitting schemes.On the finite differences grid,Δt corresponds to the time step, and hx,hy are the grid sizes in the different spatial
directions. The time nΔt is denoted by tn,
and i,j refer to the spatial coordinates of the grid
point (ihx,jhy).
Let un denote the grid function on the time level n,
and ui,jn be the specific value of un at point i,j.In Section3.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 pointsxj=jh,yk=kh,tn=nΔt, where h>0 is the spatial grid size and Δt the time step. Defining the grid function Uj,kn=U(xj,yk,tn),
the basic explicit scheme is(3.5)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ2Uj,kn+fj,kn,where ℳ2 is a difference operator(3.6)ℳ2=((λ+2μ)Dx2+μDy2(λ+μ)D0xD0y(λ+μ)D0xD0y(λ+2μ)Dy2+μDx2),
and we use the standard
difference operator notation:
(3.7)D+xvj,k=1h(vj+1,k−vj,k),D−xvj,k=D+xvj−1,k,D0x=12(D+x+D−x),Dx2=D+xD−xℳ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](3.8)Δt<hλ+3μ.Replacingℳ2 with ℳ4,
a fourth-order difference operator given by(3.9)ℳ4=((λ+2μ)(1−h212Dx2)Dx2+μ(1−h212Dy2)Dy2(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+2μ)(1−h212Dy2)Dy2+μ(1−h212Dx2)Dx2),
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
(3.10)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ4Uj,kn+fj,kn+Δt212(ℳ22Uj,kn+ℳ2fi,jn+∂ttfi,jn),where ℳ22 is a second-order approximation to the squared
right-hand side operator in (2.4a). As it only needs to be second-order accurate, ℳ22 has the same extent in space as ℳ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:(3.11)ρUj,kn+1−2Uj,kn+Uj,kn−1Δt2=ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)+θfj,kn+1+(1−2θ)fj,kn+θfj,kn−1.
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(4.1)∂ttu(t)=(A+B+C)u(t)+f(t),t∈(tn,tn+1),u(tn)=u0,u′(tn)=u1,
where the initial functions u0 and u1 are given. We could also apply for u1 that u′(tn)=(u(tn)−u(tn−1))/Δt+O(Δt)=u1.
Consequently, we have u(tn−1)≈u0−Δtu1.
The right-hand side f(t) is given as a force term.The spatial discretization terms are given
by(4.2)A=∂2∂x2,B=∂2∂y2,C=∂2∂z2,
where the approximated
discretization is given by(4.3)Au(x,y,z)≈u(x+Δx,y,z)−2u(x,y,z)+u(x−Δx,y,z)Δx2,Bu(x,y,z)≈u(x,y+Δy,z)−2u(x,y,z)+u(x,y−Δy,z)Δy2,Cu(x,y,z)≈u(x,y,z+Δz)−2u(x,y,z)+u(x,y,z−Δz)Δz2.We could decouple the equation into 3 simpler
equations obtaining a method of second order:(4.4a)u˜˜−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2Bu(tn)+Δt2Cu(tn)+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),(4.4b)u˜−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2B(ηu˜+(1−2η)u(tn)+ηu(tn−1))+Δt2Cu(tn)+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),(4.4c)u(tn+1)−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2B(ηu˜+(1−2η)u(tn)+ηu(tn−1))+Δt2C(ηu(tn+1)+(1−2η)u(tn)+ηu(tn−1))+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),  where the result is given as u(tn+1) with the initial conditions u(tn)=u0,u(tn−1)=u0−Δtu1, 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(tn+1).The underlying idea consists of the approximation of
the pairwise operators:(4.5)Δt2Aη(u˜˜−2u(tn)+u(tn−1))≈0,Δt2Bη(u˜−2u(tn)+u(tn−1))≈0,
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(4.6)∂ttu(t)=(A+B+C)u(t)+h(t),t∈(tn,tn+1),u(tn)=g(t),u′(tn)=f(t),
where A=D1(x,y,z)(∂2/∂x2),B=D2(x,y,z)(∂2/∂y2),
and C=D3(x,y,z)(∂2/∂z2) are the spatial operators. The
wave-propagation functions are as follows: (4.7)D1(x,y,z),D2(x,y,z),D3(x,y,z):ℝ3→ℝ+.Hence, for 3 operators, we have the following
second-order splitting method:(4.8)u˜˜−2u˜˜(tn)+u˜˜(tn−1)=Δt2A(ηu˜˜+(1−2η)u˜˜(tn)+ηu˜˜(tn−1))+Δt2Bu˜˜(tn)+Δt2Cu˜˜(tn)+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)),u˜−2u˜(tn)+u˜(tn−1)=Δt2A(ηu˜˜+(1−2η)u˜(tn)+ηu˜(tn−1))+Δt2B(ηu˜+(1−2η)u˜(tn)+ηu˜(tn−1))+Δt2Cu˜(tn)+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)),u(tn+1)−2u^(tn)+u^(tn−1)=Δt2A(ηu˜˜+(1−2η)u^(tn)+ηu^(tn−1))+Δt2B(ηu˜+(1−2η)u^(tn)+ηu^(tn−1))+Δt2C(ηu(tn+1)+(1−2η)u^(tn)+ηu^(tn−1))+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)), where the result is given as u(tn+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:(4.9)∂u˜˜∂n=0
is split in(4.10)∂u˜˜∂xnx+∂u˜˜∂yny+∂u˜˜∂znz=0;∂u˜∂n=0
is split in
(4.11)∂u˜˜∂xnx+∂u˜∂yny+∂u˜∂znz=0;∂u(tn+1)∂n=0
is split in
(4.12)∂u˜˜∂xnx+∂u˜∂yny+∂un+1∂znz=0.(3) Outflowing values, we have to decouple the values
into the different directions:(4.13)nD∇u˜˜=uout
is split in
(4.14)D1∂xu˜˜nx+D2∂yu˜˜ny+D3∂zu˜˜nz=uout;nD∇u˜=uout
is split in
(4.15)D1∂xu˜˜nx+D2∂yu˜ny+D3∂zu˜nz=uout;nD∇un+1=uout
is split in
(4.16)D1∂xu˜˜nx+D2∂yu˜ny+D3∂zun+1nz=uout,
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)u(tn)=u0,u(tn−1)=u0−Δtu1+Δt22((A+B+C)u0+f(tn))+O(Δt3),u(tn−1)=u0−Δtu1+Δt22((A+B+C)(u0−Δt3u1+Δt212(A+B+C)u0))+Δt22f(tn)−Δt36∂f(tn)∂t+Δt424∂2f(tn)∂t2+O(Δt5).Remark 4.1.
By solving the two or three
splitting steps, it is important to mention that each solutionu˜˜,u˜,
and u 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
(4.18)un+1,0−2un+un−1=Δt2(A+B)un,un+1,1−un+1,0=Δt2ηA(un+1−2un+un−1),un+1−un+1,1=Δt2ηB(un+1−2un+un−1), where η∈(0.0,0.5) and A,B are the spatial discretized operators.If we eliminate the intermediate values in
(4.18), we obtain
(4.19)un+1−2un+un−1=Δt2(A+B)(ηun+1+(1−2η)un+ηun−1)+Bη(un+1−2un+un−1),
where Bη=η2Δt4(AB) and thus Bη(un+1−2un+un−1)=O(Δt6).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 ofO(hp),p=2,4,
for the discretization in space, where h=hx=hy 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(4.20)utt−Au−(∂tt¯u−A˜u)=O(Δt2),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]:(4.21)∂tt¯un−A˜(θun+1+(1−2θ)un+θun−1)−B˜(un+1−2un+un−1)=0,
where B˜=θ2Δt2A˜1A˜2.
Therefore, we obtain a splitting error ofB˜(un+1−2un+un−1).
Sufficient smoothness assumed
that we have(un+1−2un+un−1)=O(Δt2),
and we obtain B˜(un+1−2un+un−1)=O(Δt4).
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(4.22)∥(1−Δt2B˜)1/2∂t+un∥2+𝒫+(un,θ)≤∥(1−Δt2B˜)1/2∂t+u0∥2+𝒫+(u0,θ),
where θ∈[0.25,0.5] and 𝒫±(uj,θ):=θ(A˜uj,uj)+θ(A˜uj±1,uj±1)+(1−2θ)(A˜uj,uj±1).Proof.
We have to proof the theorem for a test function∂¯tun,
where ∂¯t denotes the central difference.
Forj≥1, we have(4.23)((1−Δt2B˜)∂¯ttuj,∂¯tuj)+(A˜(θuj+1−(1−2θ)uj+θuj−1),∂¯tuj)=0.
Multiplying withΔt and summarizing over j yield(4.24)∑j=1n((1−Δt2B˜)∂¯ttuj,∂¯tuj)Δt+(A˜((θuj+1−(1−2θ)uj+θuj−1),∂¯tuj),∂¯tuj)Δt=0.
We can derive the identities
(4.25)((1−Δt2B˜)∂¯ttuj,∂¯tuj)Δt=12∥(1−Δt2B˜)1/2∂t+uj∥2−12∥(1−Δt2B˜)1/2∂t−uh∥2,(A˜(θuj+1−(1−2θ)uj+θuj−1),∂¯tuj)Δt=12(𝒫+(uj,θ)−𝒫−(uj,θ)),
and obtain the result
(4.26)∥(1−Δt2B˜)1/2∂t+un∥2+𝒫+(un,θ)≤∥(1−Δt2B˜)1/2∂t+u0∥2+𝒫+(u0,θ),
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 multiplierA˜1A˜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].
(2) The truncation error isO(Δt2+hp),p≥2, for θ∈[0,0.5].
(3)θ=1/12, we have a fourth-order method in time O(Δt2+hp),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=Δt2/Δxmax2Dmax, where xmax are the maximal spatial step and Dmax 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 ℳ4 into three parts: ℳxx,ℳyy,
and ℳxy,
where we have
(4.27)ℳxx=((λ+2μ)(1−h212Dx2)Dx200μ(1−h212Dx2)Dx2),ℳyy=(μ(1−h212Dy2)Dy200(λ+2μ)(1−h212Dy2)Dy2),ℳxy=ℳ4−ℳxx−ℳyy.
Our proposed split method has
the following steps:
(4.28)(1)ρUj,k*−2Uj,kn+Uj,kn−1Δt2=ℳ4Uj,kn+θfj,kn+1+(1−2θ)fj,kn+θfj,kn−1,(2)ρUj,k**−Uj,k*Δt2=θℳxx(Uj,k**−2Uj,kn+Uj,kn−1)+θ2ℳxy(Uj,k*−2Uj,kn+Uj,kn−1),(3)ρUj,kn+1−Uj,k**Δt2=θℳyy(Uj,kn+1−2Uj,kn+Uj,kn−1)+θ2ℳxy(Uj,k**−2Uj,kn+Uj,kn−1).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 functionsU(x,t) the following discretization order:
(4.29)ℳ4U=μ∇2U+(λ+μ)∇(∇⋅U)+𝒪(h4).
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 solutionsU is 𝒪(Δt4),
where it is assumed that Δt=𝒪(h).Proof.
We assume in the following thatf=(0,0)T.
We add (4.28) and obtain, like in the scalar case [1], the following result for the discretized
equations
(4.30)D+tD−tUj,kn−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)−𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
where 𝒩4,θ=θ2Δt2(ℳxxℳyy+ℳxxℳxy+ℳxyℳyy)+θ3Δt4ℳxxℳyyℳxy.
We, therefore, obtain a splitting error of𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1).
For sufficient smoothness, we have(Uj,kn+1−2Uj,kn+Uj,kn−1)=𝒪(Δt2) and we obtain 𝒩4,θ(Un+1−2Un+Un−1)=𝒪(Δt4).
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 theL2(Ω).In the following, we introduce the notation of the
norms.Remark 4.8.
For our system, we extend theL2-norm as
(4.31)∥U∥L22=(U,U)L2=∫ΩU2dx,
where U2=u2+v2 or U2=u2+v2+w2 in two and three dimensions.Remark 4.9.
For a discrete grid functionUj,kn,
the L2-norm is given as(4.32)∫Ω(Ujkn)2dx=Δx2∑j,kMUj,kn,
where Δx is the uniform grid length in x and y,M is the number of grid points in x and y.
Further, Uj,kn is the solution at grid point (j,k) and at time tn.Remark 4.10.
The matrix(4.33)𝒩4,θ=θ2Δt2(ℳxxℳyy+ℳxxℳxy+ℳxyℳyy)+θ3Δt4ℳxxℳyyℳxy,
where ℳxx,ℳyy, and ℳxy 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.34)(𝒩4,θU,U)L2=∫Ω(𝒩4,θUU)dx.
Consequently, we can
denote
(4.35)∥𝒩4,θU∥≤ω∥U∥∀U∈Hd,
where ω∈ℝ+ is the weight and 𝒩4,θ is bounded. d is the dimension, and H 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
(4.36)∥(1−Δt2ω˜)1/2D+tUj,kn∥2+𝒫+(Uj,kn,θ)≤∥(1−Δt2ω˜)1/2D+tUj,k0∥2+𝒫+(Uj,k0,θ),
where one
has 𝒫±(Uj,kn,θ):=θ(ℳ4Uj,kn,Uj,kn)+θ(ℳ4Uj,kn±1,Uj,kn±1)+(1−2θ)(ℳ4Uj,kn,Uj,kn±1) and 𝒫±≥0 for θ∈[0.25,0.5].
Further, 1−Δt2ω˜∈ℝ+ is the factor for the weighted norm (ℐ−Δt2𝒩4,θ)U≤ω˜U for all U∈Hd.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-functionD0tUj,kn.
The following result is given for the discretized
equations, see also (4.30):(4.37)(ℐ−Δt2𝒩4,θ)D+tD−tUj,kn−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)=0.
So forn≥1,
we can rewrite (4.37) for the stability proof as
follows:
(4.38)((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)−(ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)=0.
Multiplying withΔt and summarizing over the time levels, we
obtain
(4.39)∑n((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)Δt−∑n(ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)Δt=0,
for each term of the sum, one
can derive the following identities. So for ℐ−Δt2𝒩4,θ,
we have
(4.40)((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)Δt=12((ℐ−Δt2𝒩4,θ)(D+t−D−t)Uj,kn,(D+t+D−t)Uj,kn)=∫Ω((ℐ−Δt2𝒩4,θ)(D+t−D−t))T(D+t+D−t)Uj,kndx≤(1−Δt2ω˜)∫Ω(D+tUj,kn)2(D−tUj,kn)2dx, where the operator ℐ−Δt2𝒩4,θ is symmetric and positive-definite and we can
apply the weighted norm, see Remark 4.11 and [11].
We obtain the following result:(4.41)(1−Δt2ω˜)∫Ω(D+tUj,kn)2(D−tUj,kn)2dx=12∥(1−Δt2ω˜)1/2D+tUj,kn∥2−12∥(1−Δt2ω˜)1/2D−tUj,kn∥2.
Further, for−ℳ4, we have
(4.42)(−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)Δt=12(𝒫+(Uj,kn,θ)−𝒫-(Uj,kn,θ)).
Due to the result of the operators:(4.43)𝒫-(Uj,kn,θ)=𝒫+(Uj,kn−1,θ),D−tUj,kn=D+tUj,kn−1, we can recursively derive the following
result:
(4.44)∥(1−Δt2ω˜)1/2D+tUj,kn∥2+𝒫+(Uj,kn,θ)≤∥(1−Δt2ω˜)1/2D+tUj,k0∥2+𝒫+(Uj,k0,θ),
where for θ∈[0.25,0.5],
we have 𝒫+(Uj,kn,θ)≥0 for all n∈ℕ+,
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
(4.45)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)+ℳ4Uj,kn+𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
where ℳ4 is a fourth-order discretization scheme in
space.The proof is given as follows.Proof.
We consider the following Taylor-expansion:(4.46)∂ttUj,kn=D+tD−tUj,kn−Δt212∂ttttUj,kn+𝒪(Δt4).
Furthermore, we have
(4.47)∂ttttUj,kn≈ℳ4∂ttUj,kn,
and we can rewrite (4.46) as
(4.48)∂ttUj,kn≈D+tD−tUj,kn−Δt212ℳ4∂ttUj,kn+𝒪(Δt4)≈D+tD−tUj,kn−Δt212ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)+𝒪(Δt4).
So the fourth-order time-stepping algorithm can be
formulated as(4.49)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)−ℳ4Uj,kn=0.
The split method (4.28) becomes
(4.50)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)−ℳ4Uj,kn−𝒩4,1/12(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
and we obtain a fourth-order
split scheme (cf. the scalar case [1]).Remark 4.15.
As follows form Theorem4.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.


## 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(4.1)∂ttu(t)=(A+B+C)u(t)+f(t),t∈(tn,tn+1),u(tn)=u0,u′(tn)=u1,
where the initial functions u0 and u1 are given. We could also apply for u1 that u′(tn)=(u(tn)−u(tn−1))/Δt+O(Δt)=u1.
Consequently, we have u(tn−1)≈u0−Δtu1.
The right-hand side f(t) is given as a force term.The spatial discretization terms are given
by(4.2)A=∂2∂x2,B=∂2∂y2,C=∂2∂z2,
where the approximated
discretization is given by(4.3)Au(x,y,z)≈u(x+Δx,y,z)−2u(x,y,z)+u(x−Δx,y,z)Δx2,Bu(x,y,z)≈u(x,y+Δy,z)−2u(x,y,z)+u(x,y−Δy,z)Δy2,Cu(x,y,z)≈u(x,y,z+Δz)−2u(x,y,z)+u(x,y,z−Δz)Δz2.We could decouple the equation into 3 simpler
equations obtaining a method of second order:(4.4a)u˜˜−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2Bu(tn)+Δt2Cu(tn)+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),(4.4b)u˜−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2B(ηu˜+(1−2η)u(tn)+ηu(tn−1))+Δt2Cu(tn)+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),(4.4c)u(tn+1)−2u(tn)+u(tn−1)=Δt2A(ηu˜˜+(1−2η)u(tn)+ηu(tn−1))+Δt2B(ηu˜+(1−2η)u(tn)+ηu(tn−1))+Δt2C(ηu(tn+1)+(1−2η)u(tn)+ηu(tn−1))+Δt2(ηf(tn+1)+(1−2η)f(tn)+ηf(tn−1)),  where the result is given as u(tn+1) with the initial conditions u(tn)=u0,u(tn−1)=u0−Δtu1, 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(tn+1).The underlying idea consists of the approximation of
the pairwise operators:(4.5)Δt2Aη(u˜˜−2u(tn)+u(tn−1))≈0,Δt2Bη(u˜−2u(tn)+u(tn−1))≈0,
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(4.6)∂ttu(t)=(A+B+C)u(t)+h(t),t∈(tn,tn+1),u(tn)=g(t),u′(tn)=f(t),
where A=D1(x,y,z)(∂2/∂x2),B=D2(x,y,z)(∂2/∂y2),
and C=D3(x,y,z)(∂2/∂z2) are the spatial operators. The
wave-propagation functions are as follows: (4.7)D1(x,y,z),D2(x,y,z),D3(x,y,z):ℝ3→ℝ+.Hence, for 3 operators, we have the following
second-order splitting method:(4.8)u˜˜−2u˜˜(tn)+u˜˜(tn−1)=Δt2A(ηu˜˜+(1−2η)u˜˜(tn)+ηu˜˜(tn−1))+Δt2Bu˜˜(tn)+Δt2Cu˜˜(tn)+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)),u˜−2u˜(tn)+u˜(tn−1)=Δt2A(ηu˜˜+(1−2η)u˜(tn)+ηu˜(tn−1))+Δt2B(ηu˜+(1−2η)u˜(tn)+ηu˜(tn−1))+Δt2Cu˜(tn)+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)),u(tn+1)−2u^(tn)+u^(tn−1)=Δt2A(ηu˜˜+(1−2η)u^(tn)+ηu^(tn−1))+Δt2B(ηu˜+(1−2η)u^(tn)+ηu^(tn−1))+Δt2C(ηu(tn+1)+(1−2η)u^(tn)+ηu^(tn−1))+Δt2(ηh(tn+1)+(1−2η)h(tn)+ηh(tn−1)), where the result is given as u(tn+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:(4.9)∂u˜˜∂n=0
is split in(4.10)∂u˜˜∂xnx+∂u˜˜∂yny+∂u˜˜∂znz=0;∂u˜∂n=0
is split in
(4.11)∂u˜˜∂xnx+∂u˜∂yny+∂u˜∂znz=0;∂u(tn+1)∂n=0
is split in
(4.12)∂u˜˜∂xnx+∂u˜∂yny+∂un+1∂znz=0.(3) Outflowing values, we have to decouple the values
into the different directions:(4.13)nD∇u˜˜=uout
is split in
(4.14)D1∂xu˜˜nx+D2∂yu˜˜ny+D3∂zu˜˜nz=uout;nD∇u˜=uout
is split in
(4.15)D1∂xu˜˜nx+D2∂yu˜ny+D3∂zu˜nz=uout;nD∇un+1=uout
is split in
(4.16)D1∂xu˜˜nx+D2∂yu˜ny+D3∂zun+1nz=uout,
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)u(tn)=u0,u(tn−1)=u0−Δtu1+Δt22((A+B+C)u0+f(tn))+O(Δt3),u(tn−1)=u0−Δtu1+Δt22((A+B+C)(u0−Δt3u1+Δt212(A+B+C)u0))+Δt22f(tn)−Δt36∂f(tn)∂t+Δt424∂2f(tn)∂t2+O(Δt5).Remark 4.1.
By solving the two or three
splitting steps, it is important to mention that each solutionu˜˜,u˜,
and u 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
(4.18)un+1,0−2un+un−1=Δt2(A+B)un,un+1,1−un+1,0=Δt2ηA(un+1−2un+un−1),un+1−un+1,1=Δt2ηB(un+1−2un+un−1), where η∈(0.0,0.5) and A,B are the spatial discretized operators.If we eliminate the intermediate values in
(4.18), we obtain
(4.19)un+1−2un+un−1=Δt2(A+B)(ηun+1+(1−2η)un+ηun−1)+Bη(un+1−2un+un−1),
where Bη=η2Δt4(AB) and thus Bη(un+1−2un+un−1)=O(Δt6).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 ofO(hp),p=2,4,
for the discretization in space, where h=hx=hy 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(4.20)utt−Au−(∂tt¯u−A˜u)=O(Δt2),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]:(4.21)∂tt¯un−A˜(θun+1+(1−2θ)un+θun−1)−B˜(un+1−2un+un−1)=0,
where B˜=θ2Δt2A˜1A˜2.
Therefore, we obtain a splitting error ofB˜(un+1−2un+un−1).
Sufficient smoothness assumed
that we have(un+1−2un+un−1)=O(Δt2),
and we obtain B˜(un+1−2un+un−1)=O(Δt4).
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(4.22)∥(1−Δt2B˜)1/2∂t+un∥2+𝒫+(un,θ)≤∥(1−Δt2B˜)1/2∂t+u0∥2+𝒫+(u0,θ),
where θ∈[0.25,0.5] and 𝒫±(uj,θ):=θ(A˜uj,uj)+θ(A˜uj±1,uj±1)+(1−2θ)(A˜uj,uj±1).Proof.
We have to proof the theorem for a test function∂¯tun,
where ∂¯t denotes the central difference.
Forj≥1, we have(4.23)((1−Δt2B˜)∂¯ttuj,∂¯tuj)+(A˜(θuj+1−(1−2θ)uj+θuj−1),∂¯tuj)=0.
Multiplying withΔt and summarizing over j yield(4.24)∑j=1n((1−Δt2B˜)∂¯ttuj,∂¯tuj)Δt+(A˜((θuj+1−(1−2θ)uj+θuj−1),∂¯tuj),∂¯tuj)Δt=0.
We can derive the identities
(4.25)((1−Δt2B˜)∂¯ttuj,∂¯tuj)Δt=12∥(1−Δt2B˜)1/2∂t+uj∥2−12∥(1−Δt2B˜)1/2∂t−uh∥2,(A˜(θuj+1−(1−2θ)uj+θuj−1),∂¯tuj)Δt=12(𝒫+(uj,θ)−𝒫−(uj,θ)),
and obtain the result
(4.26)∥(1−Δt2B˜)1/2∂t+un∥2+𝒫+(un,θ)≤∥(1−Δt2B˜)1/2∂t+u0∥2+𝒫+(u0,θ),
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 multiplierA˜1A˜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].
(2) The truncation error isO(Δt2+hp),p≥2, for θ∈[0,0.5].
(3)θ=1/12, we have a fourth-order method in time O(Δt2+hp),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=Δt2/Δxmax2Dmax, where xmax are the maximal spatial step and Dmax 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 ℳ4 into three parts: ℳxx,ℳyy,
and ℳxy,
where we have
(4.27)ℳxx=((λ+2μ)(1−h212Dx2)Dx200μ(1−h212Dx2)Dx2),ℳyy=(μ(1−h212Dy2)Dy200(λ+2μ)(1−h212Dy2)Dy2),ℳxy=ℳ4−ℳxx−ℳyy.
Our proposed split method has
the following steps:
(4.28)(1)ρUj,k*−2Uj,kn+Uj,kn−1Δt2=ℳ4Uj,kn+θfj,kn+1+(1−2θ)fj,kn+θfj,kn−1,(2)ρUj,k**−Uj,k*Δt2=θℳxx(Uj,k**−2Uj,kn+Uj,kn−1)+θ2ℳxy(Uj,k*−2Uj,kn+Uj,kn−1),(3)ρUj,kn+1−Uj,k**Δt2=θℳyy(Uj,kn+1−2Uj,kn+Uj,kn−1)+θ2ℳxy(Uj,k**−2Uj,kn+Uj,kn−1).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 functionsU(x,t) the following discretization order:
(4.29)ℳ4U=μ∇2U+(λ+μ)∇(∇⋅U)+𝒪(h4).
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 solutionsU is 𝒪(Δt4),
where it is assumed that Δt=𝒪(h).Proof.
We assume in the following thatf=(0,0)T.
We add (4.28) and obtain, like in the scalar case [1], the following result for the discretized
equations
(4.30)D+tD−tUj,kn−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)−𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
where 𝒩4,θ=θ2Δt2(ℳxxℳyy+ℳxxℳxy+ℳxyℳyy)+θ3Δt4ℳxxℳyyℳxy.
We, therefore, obtain a splitting error of𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1).
For sufficient smoothness, we have(Uj,kn+1−2Uj,kn+Uj,kn−1)=𝒪(Δt2) and we obtain 𝒩4,θ(Un+1−2Un+Un−1)=𝒪(Δt4).
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 theL2(Ω).In the following, we introduce the notation of the
norms.Remark 4.8.
For our system, we extend theL2-norm as
(4.31)∥U∥L22=(U,U)L2=∫ΩU2dx,
where U2=u2+v2 or U2=u2+v2+w2 in two and three dimensions.Remark 4.9.
For a discrete grid functionUj,kn,
the L2-norm is given as(4.32)∫Ω(Ujkn)2dx=Δx2∑j,kMUj,kn,
where Δx is the uniform grid length in x and y,M is the number of grid points in x and y.
Further, Uj,kn is the solution at grid point (j,k) and at time tn.Remark 4.10.
The matrix(4.33)𝒩4,θ=θ2Δt2(ℳxxℳyy+ℳxxℳxy+ℳxyℳyy)+θ3Δt4ℳxxℳyyℳxy,
where ℳxx,ℳyy, and ℳxy 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.34)(𝒩4,θU,U)L2=∫Ω(𝒩4,θUU)dx.
Consequently, we can
denote
(4.35)∥𝒩4,θU∥≤ω∥U∥∀U∈Hd,
where ω∈ℝ+ is the weight and 𝒩4,θ is bounded. d is the dimension, and H 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
(4.36)∥(1−Δt2ω˜)1/2D+tUj,kn∥2+𝒫+(Uj,kn,θ)≤∥(1−Δt2ω˜)1/2D+tUj,k0∥2+𝒫+(Uj,k0,θ),
where one
has 𝒫±(Uj,kn,θ):=θ(ℳ4Uj,kn,Uj,kn)+θ(ℳ4Uj,kn±1,Uj,kn±1)+(1−2θ)(ℳ4Uj,kn,Uj,kn±1) and 𝒫±≥0 for θ∈[0.25,0.5].
Further, 1−Δt2ω˜∈ℝ+ is the factor for the weighted norm (ℐ−Δt2𝒩4,θ)U≤ω˜U for all U∈Hd.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-functionD0tUj,kn.
The following result is given for the discretized
equations, see also (4.30):(4.37)(ℐ−Δt2𝒩4,θ)D+tD−tUj,kn−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1)=0.
So forn≥1,
we can rewrite (4.37) for the stability proof as
follows:
(4.38)((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)−(ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)=0.
Multiplying withΔt and summarizing over the time levels, we
obtain
(4.39)∑n((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)Δt−∑n(ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)Δt=0,
for each term of the sum, one
can derive the following identities. So for ℐ−Δt2𝒩4,θ,
we have
(4.40)((ℐ−Δt2𝒩4,θ)D+tD−tUj,kn,D0tUj,kn)Δt=12((ℐ−Δt2𝒩4,θ)(D+t−D−t)Uj,kn,(D+t+D−t)Uj,kn)=∫Ω((ℐ−Δt2𝒩4,θ)(D+t−D−t))T(D+t+D−t)Uj,kndx≤(1−Δt2ω˜)∫Ω(D+tUj,kn)2(D−tUj,kn)2dx, where the operator ℐ−Δt2𝒩4,θ is symmetric and positive-definite and we can
apply the weighted norm, see Remark 4.11 and [11].
We obtain the following result:(4.41)(1−Δt2ω˜)∫Ω(D+tUj,kn)2(D−tUj,kn)2dx=12∥(1−Δt2ω˜)1/2D+tUj,kn∥2−12∥(1−Δt2ω˜)1/2D−tUj,kn∥2.
Further, for−ℳ4, we have
(4.42)(−ℳ4(θUj,kn+1+(1−2θ)Uj,kn+θUj,kn−1),D0tUj,kn)Δt=12(𝒫+(Uj,kn,θ)−𝒫-(Uj,kn,θ)).
Due to the result of the operators:(4.43)𝒫-(Uj,kn,θ)=𝒫+(Uj,kn−1,θ),D−tUj,kn=D+tUj,kn−1, we can recursively derive the following
result:
(4.44)∥(1−Δt2ω˜)1/2D+tUj,kn∥2+𝒫+(Uj,kn,θ)≤∥(1−Δt2ω˜)1/2D+tUj,k0∥2+𝒫+(Uj,k0,θ),
where for θ∈[0.25,0.5],
we have 𝒫+(Uj,kn,θ)≥0 for all n∈ℕ+,
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
(4.45)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)+ℳ4Uj,kn+𝒩4,θ(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
where ℳ4 is a fourth-order discretization scheme in
space.The proof is given as follows.Proof.
We consider the following Taylor-expansion:(4.46)∂ttUj,kn=D+tD−tUj,kn−Δt212∂ttttUj,kn+𝒪(Δt4).
Furthermore, we have
(4.47)∂ttttUj,kn≈ℳ4∂ttUj,kn,
and we can rewrite (4.46) as
(4.48)∂ttUj,kn≈D+tD−tUj,kn−Δt212ℳ4∂ttUj,kn+𝒪(Δt4)≈D+tD−tUj,kn−Δt212ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)+𝒪(Δt4).
So the fourth-order time-stepping algorithm can be
formulated as(4.49)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)−ℳ4Uj,kn=0.
The split method (4.28) becomes
(4.50)D+tD−tUj,kn−112ℳ4(Uj,kn+1−2Uj,kn+Uj,kn−1)−ℳ4Uj,kn−𝒩4,1/12(Uj,kn+1−2Uj,kn+Uj,kn−1)=0,
and we obtain a fourth-order
split scheme (cf. the scalar case [1]).Remark 4.15.
As follows form Theorem4.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:(5.1)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,u(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt),on∂Ω×(0,T),
where the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given by(5.2)uana(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).For the approximation error, we choose theL1-norm.TheL1-norm is given by
(5.3)errL1:=∑i,j=1,…,mVi,j|u(xi,yj,tn)−uana(xi,yj,tn)|,
where u(xi,yj,tn) is the numerical and uana(xi,yj,tn) is the analytical solution and Vi,j=ΔxΔy.Our test examples are organized as follows.(1)
The non-stiff case. We chooseD1=D2=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 unum and uana are taken from the center of our domain.(2)
The stiff case. We chooseD1=D2=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 unum and uana 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
Figure1 and their results in Figure 2.Figure 1
Numerical
error for standard and modified methods, with respect to theη parameter and Dirichlet boundary conditions.Numerical resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Dirichlet boundary conditions (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)
(classical method).
(a)(b)The numerical errors for the stiff case with Dirichlet
boundary conditions are presented in Figure3
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 (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01).Numerical
approximation and error function for the Dirichlet boundary in the stiff case (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01.).
(a)(b)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:(5.4)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,∂u(x,y,t)∂n=∂uanaly(x,y,t)∂n=0,on∂Ω×(0,T),
where Ω=[0,1]×[0,1].D1=1,D2=0.5 and the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given as(5.5)canaly(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).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 Figure5 and their results
in Figure 6.Figure 5
Numerical error
for standard and modified methods, with respect to theη parameter and Neumann boundary conditions.Numerical
resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Neumann boundary condition (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)  (classical method).
(a)(b)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(5.6)∂ttu=D1(x,y)∂xxu+D2(x,y)∂yyu,inΩ×(0,T),u(x,y,tn)=u0,∂tu(x,y,tn)=u1,on∂Ω×(0,T),u(x,y,t)=u2,on∂Ω×(0,T),
where D1(x,y)=0.1x+0.01y+0.01,D2(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.Remark 5.3.
We have assumed the following CFL condition:(5.7)Δt<0.5min(Δx,Δy)maxx,y∈Ω(D1(x,y),D2(x,y)).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Δx=1/256,Δy=1/256, and the time step Δt=1/256<0.390625.The model domain is given by a rectangle withΔx=1/16 and Δy=1/32.
The time steps are given by Δt=1/16 and 0≤η≤1.The numerical errors for the spatial-dependent
parameters with Dirichlet boundary conditions are presented in Figure7
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.Dirichlet
boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)The numerical errors for the spatial-dependent
parameters with Neumann boundary conditions are presented in Figure9 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.Neumann boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)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 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.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 timet=0,
we know the value at level n=0 as U0=g0. The value at level n=−1 can be obtained by Taylor expansion as
(5.8)U−1=U0−Δt∂tU0+Δt22∂ttU0−Δt36∂tttU0+Δt424∂ttttU0+O(Δt5),
where we use
(5.9a)∂tUj,k0=g1j,k,(5.9b)∂ttUj,k0≈1ρ(ℳ4g0j,k)+fj,k),(5.9c)∂tttUj,k0≈1ρ(ℳ4g1j,k)+∂tfj,k0),(5.9d)∂ttttUj,k0≈1ρ(ℳ22g0j,k)+ℳ4fj,k0+∂ttfj,k0), and also for (5.9c) and
(5.9d),
(5.9e)∂tfj,k0≈fj,k1−fj,k−12Δt,(5.9f)∂ttfj,k0≈fj,k1−2fj,k0+fj,k−1Δt2.
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 Un+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(5.10)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)−2Vs2sin(x)sin(t−y)−(λ+μ)(cos(t−x)cos(y)+sin(y)sin(t−y)))T,
giving the analytical solution
(5.11)Utrue=(sin(x−t)sin(y),sin(y−t)sin(x))T.
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 ∥Uj,k∥=max(maxj,k|uj,k|,maxj,k|vj,k|). As can be seen
we get the expected 4th order convergence for problems with smooth solutions.Table 1
Errors in max-norm for decreasingh and smooth analytical solution Utrue.
Convergence rate indicates a fourth-order convergence for the split
scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.051.7683e-07—2.5403e-07—0.0251.2220e-083.8552.1104e-083.5890.01257.9018e-103.9511.4376e-093.8760.0061255.0013e-113.9829.2727e-113.955To 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 decreasingh and smooth analytical solution Utrue and using the non-split scheme. Comparing with
Table 1, we see that the splitting error is very small for this
case.
t=2,eh=∥Un−Utrue∥∞hcase 1case 20.051.6878e-072.4593e-070.0251.1561e-082.0682e-080.01257.4757e-101.4205e-090.0061254.8112e-119.2573e-11


#### 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)f=Fδ(x)g(t),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
(δh2) and one 4th (δh4) order regularized approximations of the one
dimensional Dirac function,
(5.13)δh2(x˜)=1h{1+x˜,−h≤x˜<0,1−x˜,0≤x˜<h,0,elsewhere,(5.14)δh4(x˜)=1h{1+116x˜+58x˜2+16x˜3,−2h≤x˜<−h,1+12x˜−x˜2−12x˜3,−h≤x˜<0,1−12x˜−x˜2+12x˜3,0≤x˜<h,1−116x˜+x˜2−16x˜3,h≤x˜<2h,0,elsewhere,where in the above x˜=x/h.
The two and three dimensional Dirac functions are then approximated as δh2,4(x˜)δh2,4(y˜) and δh2,4(x˜)δh2,4(y˜)δh2,4(z˜).
The chosen time dependence was a smooth function given by
(5.15)g(t)={exp(−1t(1−t)),0≤t<1,0,elsewhere,
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 δh2 and δh4.
As can be seen the error is smooth and converges a small distance away from the
source. However, using δh2 limits the convergence to 2nd order, while
using δh4 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 δh2 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 t=1.Table 3
Errors in max-norm for decreasingh and analytical solution Utrue.
Convergence rate approaches 4th order after the singular forcing term goes to
zero.
t=1.1,eh=∥Un−Utrue∥∞hlog2(e2h/eh)0.051.1788e-04—0.0251.4146e-053.05880.01251.3554e-063.38360.006251.0718e-073.66060.0031257.1890e-093.8981Figure 11
The 2-logarithm of the error along a line going
through the source point for a point force located atx=0,y=0,
and approximated in space by (5.14).
Note that the error decays as 𝒪(h4) away from the source, but not near it. The
grid sizes were h=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 atx=0,y=0,
and approximated in space by (5.13).
Note that the error only decays as 𝒪(h2) away from the source. The grid sizes were h=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(5.16)ℳ4=((λ+2μ)(1−h212Dx2)Dx2+μ(1−h212Dy2Dy2+1−h212Dz2)Dz2(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+2μ)(1−h212Dy2)Dy2+μ(1−h212Dx2)Dx2+μ(1−h212Dz2)Dz2(λ+μ)(1−h26Dz2)D0z(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dy2)D0y(1−h26Dz2)D0z(λ+2μ)(1−h212Dz2Dz2+μ(1−h212Dx2)Dx2+1−h212Dy2)Dy2),
operating on grid functions Uj,k,ln defined at grid points xj,yk,zl,tn similarly to the two dimensional case. We can
split ℳ4 into six parts; ℳxx,ℳyy,ℳzz containing the three second order directional
difference operators, and ℳxy,ℳyz,ℳxz 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:(5.17)(1)ρUj,k,l*−2Uj,k,ln+Uj,k,ln−1Δt2=ℳ4Uj,k,ln+θfj,k,ln+1+(1−2θ)fj,k,ln+θfj,k,ln−1(2)ρUj,k,l**−Uj,k,l*Δt2=θℳxx(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳxz)(Uj,k,l*−2Uj,k,ln+Uj,k,ln−1)(3)ρUj,k,l***−Uj,k,l**Δt2=θℳxx(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳyz)(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)(4)ρUj,k,ln+1−Uj,k,l***Δt2=θℳxx(Uj,k,ln+1−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxz+ℳyz)(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1).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(5.18)f=−(−1+λ+4μ)sin(t−x)sin(y)sin(z)−(λ+μ)cos(x)(2sin(t)sin(y)sin(z)+cos(t)sin(y+z)),−(−1+λ+4μ)sin(x)sin(t−y)sin(z)−(λ+μ)cos(y)(2sin(t)sin(x)sin(z)+cos(t)sin(x+z)),−(λ+μ)cos(t−y)cos(z)sin(x)−sin(y)((λ+μ)cos(t−x)cos(z)+(−1+λ+4μ)sin(x)sin(t−z)))T,
giving the analytical solution
(5.19)Utrue=(sin(x−t)sin(y)sin(z),sin(y−t)sin(x)sin(z),sin(z−t)sin(x)sin(y))T.
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 t=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 decreasingh and smooth analytical solution Utrue.
Convergence rate indicates 4th order convergence for the three dimensional
split scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.14.2986e-07—1.8542e-06—0.053.5215e-083.611.3605e-073.770.0253.0489e-093.538.0969e-094.070.01252.0428e-103.904.7053e-104.10


## 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:(5.1)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,u(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt),on∂Ω×(0,T),
where the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given by(5.2)uana(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).For the approximation error, we choose theL1-norm.TheL1-norm is given by
(5.3)errL1:=∑i,j=1,…,mVi,j|u(xi,yj,tn)−uana(xi,yj,tn)|,
where u(xi,yj,tn) is the numerical and uana(xi,yj,tn) is the analytical solution and Vi,j=ΔxΔy.Our test examples are organized as follows.(1)
The non-stiff case. We chooseD1=D2=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 unum and uana are taken from the center of our domain.(2)
The stiff case. We chooseD1=D2=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 unum and uana 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
Figure1 and their results in Figure 2.Figure 1
Numerical
error for standard and modified methods, with respect to theη parameter and Dirichlet boundary conditions.Numerical resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Dirichlet boundary conditions (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)
(classical method).
(a)(b)The numerical errors for the stiff case with Dirichlet
boundary conditions are presented in Figure3
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 (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01).Numerical
approximation and error function for the Dirichlet boundary in the stiff case (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01.).
(a)(b)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:(5.4)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,∂u(x,y,t)∂n=∂uanaly(x,y,t)∂n=0,on∂Ω×(0,T),
where Ω=[0,1]×[0,1].D1=1,D2=0.5 and the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given as(5.5)canaly(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).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 Figure5 and their results
in Figure 6.Figure 5
Numerical error
for standard and modified methods, with respect to theη parameter and Neumann boundary conditions.Numerical
resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Neumann boundary condition (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)  (classical method).
(a)(b)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(5.6)∂ttu=D1(x,y)∂xxu+D2(x,y)∂yyu,inΩ×(0,T),u(x,y,tn)=u0,∂tu(x,y,tn)=u1,on∂Ω×(0,T),u(x,y,t)=u2,on∂Ω×(0,T),
where D1(x,y)=0.1x+0.01y+0.01,D2(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.Remark 5.3.
We have assumed the following CFL condition:(5.7)Δt<0.5min(Δx,Δy)maxx,y∈Ω(D1(x,y),D2(x,y)).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Δx=1/256,Δy=1/256, and the time step Δt=1/256<0.390625.The model domain is given by a rectangle withΔx=1/16 and Δy=1/32.
The time steps are given by Δt=1/16 and 0≤η≤1.The numerical errors for the spatial-dependent
parameters with Dirichlet boundary conditions are presented in Figure7
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.Dirichlet
boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)The numerical errors for the spatial-dependent
parameters with Neumann boundary conditions are presented in Figure9 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.Neumann boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)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.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:(5.1)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,u(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt),on∂Ω×(0,T),
where the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given by(5.2)uana(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).For the approximation error, we choose theL1-norm.TheL1-norm is given by
(5.3)errL1:=∑i,j=1,…,mVi,j|u(xi,yj,tn)−uana(xi,yj,tn)|,
where u(xi,yj,tn) is the numerical and uana(xi,yj,tn) is the analytical solution and Vi,j=ΔxΔy.Our test examples are organized as follows.(1)
The non-stiff case. We chooseD1=D2=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 unum and uana are taken from the center of our domain.(2)
The stiff case. We chooseD1=D2=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 unum and uana 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
Figure1 and their results in Figure 2.Figure 1
Numerical
error for standard and modified methods, with respect to theη parameter and Dirichlet boundary conditions.Numerical resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Dirichlet boundary conditions (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)
(classical method).
(a)(b)The numerical errors for the stiff case with Dirichlet
boundary conditions are presented in Figure3
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 (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01).Numerical
approximation and error function for the Dirichlet boundary in the stiff case (Δx=Δy=1/32,Δt=1/64,D1=1,D2=0.01.).
(a)(b)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:(5.4)∂ttu=D12∂xxu+D22∂yyu,inΩ×(0,T),u(x,y,0)=u0(x,y)=sin(1D1πx)sin(1D2πy),onΩ,∂tu(x,y,0)=u1(x,y)=0,onΩ,∂u(x,y,t)∂n=∂uanaly(x,y,t)∂n=0,on∂Ω×(0,T),
where Ω=[0,1]×[0,1].D1=1,D2=0.5 and the initial conditions can be written as u(x,y,tn)=u0(x,y) and u(x,y,tn−1)=u(x,y,tn+1)=u(x,y,Δt).The analytical solution is given as(5.5)canaly(x,y,t)=sin(1D1πx)sin(1D2πy)cos(2πt).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 Figure5 and their results
in Figure 6.Figure 5
Numerical error
for standard and modified methods, with respect to theη parameter and Neumann boundary conditions.Numerical
resolution of the wave equation: numerical approximation (a) and error
functions (b) for the Neumann boundary condition (Δx=Δy=1/32,Δt=1/64,D1=1,D2=1)  (classical method).
(a)(b)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(5.6)∂ttu=D1(x,y)∂xxu+D2(x,y)∂yyu,inΩ×(0,T),u(x,y,tn)=u0,∂tu(x,y,tn)=u1,on∂Ω×(0,T),u(x,y,t)=u2,on∂Ω×(0,T),
where D1(x,y)=0.1x+0.01y+0.01,D2(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.Remark 5.3.
We have assumed the following CFL condition:(5.7)Δt<0.5min(Δx,Δy)maxx,y∈Ω(D1(x,y),D2(x,y)).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Δx=1/256,Δy=1/256, and the time step Δt=1/256<0.390625.The model domain is given by a rectangle withΔx=1/16 and Δy=1/32.
The time steps are given by Δt=1/16 and 0≤η≤1.The numerical errors for the spatial-dependent
parameters with Dirichlet boundary conditions are presented in Figure7
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.Dirichlet
boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)The numerical errors for the spatial-dependent
parameters with Neumann boundary conditions are presented in Figure9 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.Neumann boundary condition: numerical solution and error function for the
spatial-dependent test example.
(a)(b)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 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.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 timet=0,
we know the value at level n=0 as U0=g0. The value at level n=−1 can be obtained by Taylor expansion as
(5.8)U−1=U0−Δt∂tU0+Δt22∂ttU0−Δt36∂tttU0+Δt424∂ttttU0+O(Δt5),
where we use
(5.9a)∂tUj,k0=g1j,k,(5.9b)∂ttUj,k0≈1ρ(ℳ4g0j,k)+fj,k),(5.9c)∂tttUj,k0≈1ρ(ℳ4g1j,k)+∂tfj,k0),(5.9d)∂ttttUj,k0≈1ρ(ℳ22g0j,k)+ℳ4fj,k0+∂ttfj,k0), and also for (5.9c) and
(5.9d),
(5.9e)∂tfj,k0≈fj,k1−fj,k−12Δt,(5.9f)∂ttfj,k0≈fj,k1−2fj,k0+fj,k−1Δt2.
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 Un+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(5.10)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)−2Vs2sin(x)sin(t−y)−(λ+μ)(cos(t−x)cos(y)+sin(y)sin(t−y)))T,
giving the analytical solution
(5.11)Utrue=(sin(x−t)sin(y),sin(y−t)sin(x))T.
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 ∥Uj,k∥=max(maxj,k|uj,k|,maxj,k|vj,k|). As can be seen
we get the expected 4th order convergence for problems with smooth solutions.Table 1
Errors in max-norm for decreasingh and smooth analytical solution Utrue.
Convergence rate indicates a fourth-order convergence for the split
scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.051.7683e-07—2.5403e-07—0.0251.2220e-083.8552.1104e-083.5890.01257.9018e-103.9511.4376e-093.8760.0061255.0013e-113.9829.2727e-113.955To 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 decreasingh and smooth analytical solution Utrue and using the non-split scheme. Comparing with
Table 1, we see that the splitting error is very small for this
case.
t=2,eh=∥Un−Utrue∥∞hcase 1case 20.051.6878e-072.4593e-070.0251.1561e-082.0682e-080.01257.4757e-101.4205e-090.0061254.8112e-119.2573e-11


### 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)f=Fδ(x)g(t),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
(δh2) and one 4th (δh4) order regularized approximations of the one
dimensional Dirac function,
(5.13)δh2(x˜)=1h{1+x˜,−h≤x˜<0,1−x˜,0≤x˜<h,0,elsewhere,(5.14)δh4(x˜)=1h{1+116x˜+58x˜2+16x˜3,−2h≤x˜<−h,1+12x˜−x˜2−12x˜3,−h≤x˜<0,1−12x˜−x˜2+12x˜3,0≤x˜<h,1−116x˜+x˜2−16x˜3,h≤x˜<2h,0,elsewhere,where in the above x˜=x/h.
The two and three dimensional Dirac functions are then approximated as δh2,4(x˜)δh2,4(y˜) and δh2,4(x˜)δh2,4(y˜)δh2,4(z˜).
The chosen time dependence was a smooth function given by
(5.15)g(t)={exp(−1t(1−t)),0≤t<1,0,elsewhere,
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 δh2 and δh4.
As can be seen the error is smooth and converges a small distance away from the
source. However, using δh2 limits the convergence to 2nd order, while
using δh4 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 δh2 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 t=1.Table 3
Errors in max-norm for decreasingh and analytical solution Utrue.
Convergence rate approaches 4th order after the singular forcing term goes to
zero.
t=1.1,eh=∥Un−Utrue∥∞hlog2(e2h/eh)0.051.1788e-04—0.0251.4146e-053.05880.01251.3554e-063.38360.006251.0718e-073.66060.0031257.1890e-093.8981Figure 11
The 2-logarithm of the error along a line going
through the source point for a point force located atx=0,y=0,
and approximated in space by (5.14).
Note that the error decays as 𝒪(h4) away from the source, but not near it. The
grid sizes were h=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 atx=0,y=0,
and approximated in space by (5.13).
Note that the error only decays as 𝒪(h2) away from the source. The grid sizes were h=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.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 timet=0,
we know the value at level n=0 as U0=g0. The value at level n=−1 can be obtained by Taylor expansion as
(5.8)U−1=U0−Δt∂tU0+Δt22∂ttU0−Δt36∂tttU0+Δt424∂ttttU0+O(Δt5),
where we use
(5.9a)∂tUj,k0=g1j,k,(5.9b)∂ttUj,k0≈1ρ(ℳ4g0j,k)+fj,k),(5.9c)∂tttUj,k0≈1ρ(ℳ4g1j,k)+∂tfj,k0),(5.9d)∂ttttUj,k0≈1ρ(ℳ22g0j,k)+ℳ4fj,k0+∂ttfj,k0), and also for (5.9c) and
(5.9d),
(5.9e)∂tfj,k0≈fj,k1−fj,k−12Δt,(5.9f)∂ttfj,k0≈fj,k1−2fj,k0+fj,k−1Δt2.
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 Un+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(5.10)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)−2Vs2sin(x)sin(t−y)−(λ+μ)(cos(t−x)cos(y)+sin(y)sin(t−y)))T,
giving the analytical solution
(5.11)Utrue=(sin(x−t)sin(y),sin(y−t)sin(x))T.
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 ∥Uj,k∥=max(maxj,k|uj,k|,maxj,k|vj,k|). As can be seen
we get the expected 4th order convergence for problems with smooth solutions.Table 1
Errors in max-norm for decreasingh and smooth analytical solution Utrue.
Convergence rate indicates a fourth-order convergence for the split
scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.051.7683e-07—2.5403e-07—0.0251.2220e-083.8552.1104e-083.5890.01257.9018e-103.9511.4376e-093.8760.0061255.0013e-113.9829.2727e-113.955To 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 decreasingh and smooth analytical solution Utrue and using the non-split scheme. Comparing with
Table 1, we see that the splitting error is very small for this
case.
t=2,eh=∥Un−Utrue∥∞hcase 1case 20.051.6878e-072.4593e-070.0251.1561e-082.0682e-080.01257.4757e-101.4205e-090.0061254.8112e-119.2573e-11


## 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)f=Fδ(x)g(t),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
(δh2) and one 4th (δh4) order regularized approximations of the one
dimensional Dirac function,
(5.13)δh2(x˜)=1h{1+x˜,−h≤x˜<0,1−x˜,0≤x˜<h,0,elsewhere,(5.14)δh4(x˜)=1h{1+116x˜+58x˜2+16x˜3,−2h≤x˜<−h,1+12x˜−x˜2−12x˜3,−h≤x˜<0,1−12x˜−x˜2+12x˜3,0≤x˜<h,1−116x˜+x˜2−16x˜3,h≤x˜<2h,0,elsewhere,where in the above x˜=x/h.
The two and three dimensional Dirac functions are then approximated as δh2,4(x˜)δh2,4(y˜) and δh2,4(x˜)δh2,4(y˜)δh2,4(z˜).
The chosen time dependence was a smooth function given by
(5.15)g(t)={exp(−1t(1−t)),0≤t<1,0,elsewhere,
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 δh2 and δh4.
As can be seen the error is smooth and converges a small distance away from the
source. However, using δh2 limits the convergence to 2nd order, while
using δh4 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 δh2 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 t=1.Table 3
Errors in max-norm for decreasingh and analytical solution Utrue.
Convergence rate approaches 4th order after the singular forcing term goes to
zero.
t=1.1,eh=∥Un−Utrue∥∞hlog2(e2h/eh)0.051.1788e-04—0.0251.4146e-053.05880.01251.3554e-063.38360.006251.0718e-073.66060.0031257.1890e-093.8981Figure 11
The 2-logarithm of the error along a line going
through the source point for a point force located atx=0,y=0,
and approximated in space by (5.14).
Note that the error decays as 𝒪(h4) away from the source, but not near it. The
grid sizes were h=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 atx=0,y=0,
and approximated in space by (5.13).
Note that the error only decays as 𝒪(h2) away from the source. The grid sizes were h=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(5.16)ℳ4=((λ+2μ)(1−h212Dx2)Dx2+μ(1−h212Dy2Dy2+1−h212Dz2)Dz2(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+2μ)(1−h212Dy2)Dy2+μ(1−h212Dx2)Dx2+μ(1−h212Dz2)Dz2(λ+μ)(1−h26Dz2)D0z(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dy2)D0y(1−h26Dz2)D0z(λ+2μ)(1−h212Dz2Dz2+μ(1−h212Dx2)Dx2+1−h212Dy2)Dy2),
operating on grid functions Uj,k,ln defined at grid points xj,yk,zl,tn similarly to the two dimensional case. We can
split ℳ4 into six parts; ℳxx,ℳyy,ℳzz containing the three second order directional
difference operators, and ℳxy,ℳyz,ℳxz 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:(5.17)(1)ρUj,k,l*−2Uj,k,ln+Uj,k,ln−1Δt2=ℳ4Uj,k,ln+θfj,k,ln+1+(1−2θ)fj,k,ln+θfj,k,ln−1(2)ρUj,k,l**−Uj,k,l*Δt2=θℳxx(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳxz)(Uj,k,l*−2Uj,k,ln+Uj,k,ln−1)(3)ρUj,k,l***−Uj,k,l**Δt2=θℳxx(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳyz)(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)(4)ρUj,k,ln+1−Uj,k,l***Δt2=θℳxx(Uj,k,ln+1−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxz+ℳyz)(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1).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(5.18)f=−(−1+λ+4μ)sin(t−x)sin(y)sin(z)−(λ+μ)cos(x)(2sin(t)sin(y)sin(z)+cos(t)sin(y+z)),−(−1+λ+4μ)sin(x)sin(t−y)sin(z)−(λ+μ)cos(y)(2sin(t)sin(x)sin(z)+cos(t)sin(x+z)),−(λ+μ)cos(t−y)cos(z)sin(x)−sin(y)((λ+μ)cos(t−x)cos(z)+(−1+λ+4μ)sin(x)sin(t−z)))T,
giving the analytical solution
(5.19)Utrue=(sin(x−t)sin(y)sin(z),sin(y−t)sin(x)sin(z),sin(z−t)sin(x)sin(y))T.
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 t=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 decreasingh and smooth analytical solution Utrue.
Convergence rate indicates 4th order convergence for the three dimensional
split scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.14.2986e-07—1.8542e-06—0.053.5215e-083.611.3605e-073.770.0253.0489e-093.538.0969e-094.070.01252.0428e-103.904.7053e-104.10


## 5.3.1. The Splitting Scheme

In three dimensions a 4th order difference
approximation of the right hand side operator becomes(5.16)ℳ4=((λ+2μ)(1−h212Dx2)Dx2+μ(1−h212Dy2Dy2+1−h212Dz2)Dz2(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dx2)D0x(1−h26Dy2)D0y(λ+2μ)(1−h212Dy2)Dy2+μ(1−h212Dx2)Dx2+μ(1−h212Dz2)Dz2(λ+μ)(1−h26Dz2)D0z(1−h26Dy2)D0y(λ+μ)(1−h26Dx2)D0x(1−h26Dz2)D0z(λ+μ)(1−h26Dy2)D0y(1−h26Dz2)D0z(λ+2μ)(1−h212Dz2Dz2+μ(1−h212Dx2)Dx2+1−h212Dy2)Dy2),
operating on grid functions Uj,k,ln defined at grid points xj,yk,zl,tn similarly to the two dimensional case. We can
split ℳ4 into six parts; ℳxx,ℳyy,ℳzz containing the three second order directional
difference operators, and ℳxy,ℳyz,ℳxz 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:(5.17)(1)ρUj,k,l*−2Uj,k,ln+Uj,k,ln−1Δt2=ℳ4Uj,k,ln+θfj,k,ln+1+(1−2θ)fj,k,ln+θfj,k,ln−1(2)ρUj,k,l**−Uj,k,l*Δt2=θℳxx(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳxz)(Uj,k,l*−2Uj,k,ln+Uj,k,ln−1)(3)ρUj,k,l***−Uj,k,l**Δt2=θℳxx(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxy+ℳyz)(Uj,k,l**−2Uj,k,ln+Uj,k,ln−1)(4)ρUj,k,ln+1−Uj,k,l***Δt2=θℳxx(Uj,k,ln+1−2Uj,k,ln+Uj,k,ln−1)+θ2(ℳxz+ℳyz)(Uj,k,l***−2Uj,k,ln+Uj,k,ln−1).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(5.18)f=−(−1+λ+4μ)sin(t−x)sin(y)sin(z)−(λ+μ)cos(x)(2sin(t)sin(y)sin(z)+cos(t)sin(y+z)),−(−1+λ+4μ)sin(x)sin(t−y)sin(z)−(λ+μ)cos(y)(2sin(t)sin(x)sin(z)+cos(t)sin(x+z)),−(λ+μ)cos(t−y)cos(z)sin(x)−sin(y)((λ+μ)cos(t−x)cos(z)+(−1+λ+4μ)sin(x)sin(t−z)))T,
giving the analytical solution
(5.19)Utrue=(sin(x−t)sin(y)sin(z),sin(y−t)sin(x)sin(z),sin(z−t)sin(x)sin(y))T.
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 t=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 decreasingh and smooth analytical solution Utrue.
Convergence rate indicates 4th order convergence for the three dimensional
split scheme.
t=2eh=∥Un−Utrue∥∞hcase 1log2(e2h/eh)case 2log2(e2h/eh)0.14.2986e-07—1.8542e-06—0.053.5215e-083.611.3605e-073.770.0253.0489e-093.538.0969e-094.070.01252.0428e-103.904.7053e-104.10


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



---

*Source: 291968-2009-02-04.xml*