See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/228656327

An Explicit Method Based on the Implicit Runge-Kutta Algorithm for Solving

Wave Equations

Article  in  Bulletin of the Seismological Society of America · November 2009

DOI: 10.1785/0120080346

CITATIONS READS
34 583

4 authors, including:

Dinghui Yang Nian Wang

Tsinghua University Institut de Physique du Globe de Paris
132 PUBLICATIONS   1,948 CITATIONS    16 PUBLICATIONS   144 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Seismic propagation in porous media with low porosity and low permeability and inversion of reservoir parameters View project

Imaging methods of lithosphere structure based on the full wave equation and its applications View project

All content following this page was uploaded by Dinghui Yang on 17 May 2014.

The user has requested enhancement of the downloaded file.

 Bulletin of the Seismological Society of America, Vol. 99, No. 6, pp. 3340–3354, December 2009, doi: 10.1785/0120080346

An Explicit Method Based on the Implicit Runge–Kutta

Algorithm for Solving Wave Equations
by Dinghui Yang, Nian Wang, Shan Chen, and Guojie Song

Abstract A new explicit differentiator series method based on the implicit Runge–
Kutta method, called the IRK-DSM in brief, is developed for solving wave equations.
To develop the new algorithm, we first transform the wave equation, usually described
as a partial differential equation (PDE), into a system of first-order ordinary differential
equations (ODEs) with respect to time t. Then we use a truncated differentiator series
method of the implicit Runge–Kutta method to solve the semidiscrete ordinary dif-
ferential equations, while the high-order spatial derivatives included in the ODEs are
approximated by the local interpolation method. We analyze the theoretical properties
of the IRK-DSM, including the stability criteria for solving the 1D and 2D acoustic-
wave equations, numerical dispersion, discretizing error, and computational efficiency
when using the IRK-DSM to model acoustic-wave fields. For comparison, we also
present the stability criteria and numerical dispersion of the so-called Lax–Wendroff
correction (LWC) methods with the fourth-order and eighth-order accuracies for the
1D case. Promising numerical results show that the IRK-DSM provides a useful tool
for large-scale practical problems because it can effectively suppress numerical dis-
persions and source-noises caused by discretizing the acoustic- and elastic-wave equa-
tions when too-coarse grids are used or the models have a large velocity contrast
between adjacent layers. Theoretical analysis and numerical modeling also demon-
strate that the IRK-DSM, through combining both the implicit Runge–Kutta scheme
with good stability condition and the approximate differentiator series method, is a
robust wave-field modeling method.

Introduction
Numerical methods of solving wave equations have high-order finite-difference method is widely used because
proven to be successful and provided useful tools in explora- of its good properties, such as high accuracy and easy im-
tion seismology. Many numerical methods of synthetic seis- plementation for computer code, but it usually suffers from
mograms have been developed and successfully applied to serious numerical dispersions when a coarse grid is used or
solve practical geophysical problems. These methods include the models have a large velocity contrast between adjacent
finite-difference (FD) methods, such as high-order compact layers (Wang et al., 2002). The finite-element method can
difference or so-called Lax–Wendroff correction (LWC) deal with irregular region and curve boundary, whereas it
schemes (e.g., Kelly et al., 1976; Dablain, 1986; Virieux, needs to solve large-scale linear algebraic equations, result-
1986; Levander, 1988; Fornberg, 1990; Lele, 1992; Igel ing in both the large storage-space requirement and the
et al., 1995; Blanch and Robertsson, 1997; Geller and Ta- costly central processing unit (CPU) time. The reflectivity
keuchi, 1998; Takeuchi and Geller, 2000; and many others); method can only be applied to media with homogeneous
the finite-element method (Yang et al., 2008); the adaptive horizontal layers. The space operators of the pseudospectral
finite-element method (Turner et al., 1956; Whiteman, method fit the Nyquist frequency, but it needs Fourier trans-
1975; Johnson, 1990; Eriksson and Johnson, 1991); the re- form that is time-consuming and the usage of a global
flectivity method (Fuchs and Müller, 1971; Booth and basis, resulting in inaccuracies in wave fields for models
Crampin, 1983a, b); the pseudospectral method (PSM) with strong heterogeneity or sharp boundaries (Komatitsch
(Kosloff and Baysal, 1982; Kosloff et al., 1984; Vlastos et al., 2000).
et al., 2003); and the spectral element method (Priolo Numerical dispersion is an important issue in numerical
and Seriani, 1991; Priolo et al., 1994; Komatitsch and Vi- seismic simulations and has been widely studied by many
lotte, 1998; Komatitsch et al., 2000), and so on. Each meth- researchers (Alford et al., 1974; Sei and Symes, 1994; Fei
od has its advantages and disadvantages. For example, the and Larner, 1995; Zhang et al., 1999; Yang et al., 2002).

3340
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3341

Roughly speaking, numerical dispersion is an unphysical ODE and the differentiator series method that was initially
phenomenon caused by discretizing the wave equation (Sei proposed to solve Cauchy problems or the linear ODE or
and Symes, 1994; Yang et al., 2002). Such an unphysical PDE (Ke, 1996; Ke and Xie, 1999). To do this, we first trans-
phenomenon affects our recognition of seismic-wave propa- form the wave equations into a system of semidiscrete ODEs,
gation. A natural idea in dealing with numerical dispersion is and then change the diagonal implicit Runge–Kutta method
to use higher accurate methods with local difference operator for solving the semidiscrete ODEs into an explicit algorithm
(Yang et al., 2006). However, the numerical dispersion or (IRK-DSM) via the truncated differentiator series. On the
undesirable ripple affects also the performance of the so- basis of such a structure, the IRK-DSM has good computa-
called high-order FD methods. For example, the tenth-order tional stability and can effectively suppress the numerical
compact FD schemes (e.g., Dablain, 1986; Wang et al., dispersion and source-noises caused by discretizing the wave
2002) that usually use more grids than low-order schemes equations when too-coarse grids are used. This allows us to
also suffer from numerical dispersions. The demand of more use larger time and spatial increments in the IRK-DSM for
grids in high-order FD methods reduces the efficiency of the modeling seismic-wave fields, resulting in reducing compu-
algorithm for parallel implementation and prevents it from tational CPU time and saving the storage for computer code.
artificial boundary treatment such as the PSM. The flux-
corrected transport (FCT) technique can lessen or eliminate
the numerical dispersion, but it is unable to fully recover lost The Explicit IRK-DSM Based on the Implicit
resolution caused by suppressing the numerical dispersion Runge–Kutta Method
when a too-coarse grid is used or the geological models have
a large velocity contrast between adjacent layers (Fei and Transform of Wave Equations
Larner, 1995; Yang et al., 2002; Zheng et al., 2006). The In a heterogeneous elastic medium, the wave equation
so-called nearly analytic discrete method (NADM) and its can be written as:
improved versions were developed by Yang et al. (2003,
2007) for reducing the numerical dispersion caused by ∂ 2U
the discretization of acoustic- and elastic-wave equations. ρ
∇ · σ  f σ
C∶ε
∂t2
These methods use the wave displacement, the velocity, and
1
their gradient fields simultaneously to restructure the wave ε
∇U  ∇UT ; (1)
displacement-fields. Hence, they can effectively suppress the 2
numerical dispersion. However, the NADM and its improved
where ρ denotes the density; U
u1 ; u2 ; u3 T , the displace-
version have a relatively small Courant number defined by
ment vector; σ and ε, the second-order symmetric stress
α
c0 Δt=Δx with the acoustic velocity c0 (Dablain, 1986;
and strain tensors, respectively; C, the fourth-order stiffness
Sei and Symes, 1994), which limits the stability range of nu-
tensor; and f
f1 ; f2 ; f3 T , the external source force.
merical calculations, whereas the optimal NADM (ONADM)
We can rewrite equation (1) as follows:
(Yang et al., 2006) cannot be applied to the two-phase porous
wave equations including the so-called particle velocity
∂U=∂t of the wave displacement U because the ONADM ∂ 2U
ρ
D · U  f; (2)
does not compute the velocity fields. ∂t2
The stability criterion α
c0 Δt=Δx ≤ αmax (αmax
denotes the maximal Courant number) gives the relationship where D is the second-order partial differential operator.
between temporal step, spatial increment, and wave velocity to For example, in the 2D anisotropic case D is defined by
keep the numerical calculation stable, and it provides the the-


oretical guide for numerical calculations. It is well-known that ∂ ∂ ∂ ∂ ∂ ∂
D
C1  C2  C3  C4 ; (3)
implicit methods, such as alternative direction implicit method ∂x ∂x ∂z ∂z ∂x ∂z
for solving PDE or diagonal implicit Runge–Kutta method
(IRK) for solving ODE, have promising stability criteria in where
theory. However, we usually need to solve large-scale linear
algebraic equations when we use the implicit methods to solve 2 3 2 3
c11 c16 c15 c15 c14 c13
wave equations, resulting in large storage requirements and 6 7 6 7
C1
4 c16 c66 c56 5; C2
4 c56 c46 c36 5;
increasing the CPU time. The explicit method has the limita-
tion of stability conditions, but it is easy to be implemented; its c15 c56 c55 c55 c45 c35
2 3 2 3
computational speed is fast at each time-advancing step. c15 c56 c55 c55 c45 c35
An explicit algorithm based on an implicit method may have 6 7 6 7
C3
4 c14 c46 c45 5; C4
4 c45 c44 c34 5; (4)
advantages of both explicit and implicit methods.
The main purpose of this article is to propose an explicit c13 c36 c35 c35 c34 c33
Runge–Kutta differentiator series method (IRK-DSM), based
on both the implicit Runge–Kutta method for solving the and cij x; z are the elastic constants.
3342 D. Yang, N. Wang, S. Chen, and G. Song

∂u1 ∂u2 ∂u3 T

Let W
∂U
∂t
 ∂t ; ∂t ; ∂t  ; then equation (2) can be K
ni;j
L

V
ni;j  1  2rΔtK ni;j  rΔtK
ni;j 
written as:
 Fi;j tn  1  rΔt; (11)
∂U ∂W 1 1 p

W
D · U  f: (5) where r
12  63.
∂t ∂t ρ ρ
Obviously, equations (10) and (11) imply that we need
Let V
U; WT ; we can further rewrite equation (5) as the to solve two systems of linear algebraic equations at each
following step of time advancing when we use the implicit algorithm
based on equations (9)–(11) to compute V
n1
i;j , resulting in an
∂V increment of the computational cost. To avoid solving the

L · V  F; (6)
∂t system of linear equations in the implicit method, we pro-
pose to use an explicit scheme in the IRK method.
where We rewrite equation (10) as follows
 

0 I 0
L
and F
:
ni;j
L
· V
ni;j  Fi;j tn  rΔt:
I  rΔtLK
ρD 0 ρf
1 1 (12)

If krΔtLk

< 1, then using the Taylor expansion, we have
From equation (6), we have the following equations:

1 L
· V
ni;j  Fi;j tn  rΔt
K ni;j
I  rΔtL
∂ 2V ∂2V X
∞

L · Vx
L · Vz; (7)
m L
· V
ni;j  Fi;j tn  rΔt:
∂t∂x ∂t∂z
rΔtL (13)
m
0

where V x
∂V=∂x and V z
∂V=∂z.

Using the truncated Taylor series to approximate the right-
Define V

V; V x ; V z T , F

F; 0; 0T , and L


hand side of equation (13), we can further obtain the follow-
diagL; L; L. From equations (6) and (7), we derive the fol-
ing equation
lowing equation:
X
2
∂ V
K ni;j

m L
· V
ni;j  Fi;j tn  rΔt:
rΔtL (14)

L
· V
 F:

(8)
∂t m
0

Obviously, calculating K ni;j in equation (14) is an explicit

computational formula.
Formulation of the Explicit IRK-DSM Similarly, we can obtain the following approximation,
deriving from equation (11),
Owing to the local elastic property of rocks, we use the
local interpolation method (Kondoh et al., 1994; Yang et al., X
2
2003, 2007) to approximate the high-order derivatives K
ni;j

V
ni;j  1  2rΔtK ni;j 

m fL
rΔtL
∂ kl U=∂xk ∂zl ni;j and ∂ kl W=∂xk ∂zl ni;j 2 ≤ k  l ≤ 3, m
0
including in the right-hand side of equation (8) through using  Fi;j tn  1  rΔtg: (15)
the values of the wave displacement, the particle velocity,
and their gradients at the grid point (i, j) and its neighbor-
ing grid points. These computational formulae of approxi- Combining equation (9) with equations (14) and (15), we
mating the second- and third-order derivatives are listed in obtain the explicit IRK-DSM. In this IRK-DSM, equations (14)
Appendix A. After discretizing the high-order spatial deriva- and (15) involve the third-order operator L
3. It means
tives at the right-hand side of equation (8), equation (8) be- that we need to compute the high-order derivatives
comes a semidiscrete ODE. ∂ kl V=∂x

k
∂zl ni;j k  l ≥ 4 as we use formulae (14) and
Now let us take a closer look at the semidiscrete ODE (15) to compute K ni;j and K
ni;j . To avoid computing these
(8). In this article, we use the diagonal implicit Runge–Kutta high-order spatial derivatives, we use a split-operator method
(IRK) method (Hairer et al., 1993) to solve the ODE (8) as to compute L
3 · V
ni;j and L
3 · V
ni;j  1  2rΔtK ni;j  . To do
follows this, we rewrite equations (14) and (15) as follows
n1 Δt n
V
i;j
V
ni;j  K i;j  K
ni;j ; (9)
2 n n
K ni;j
V

i;j  rΔtL
2 · V
ni;j  rΔt2 L
2 · V

i;j
X
2

m · Fi;j tn  rΔt;
rΔtL (16)

V
ni;j  rΔtK ni;j   Fi;j tn  rΔt;
K ni;j
L (10) m
0
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3343

X
2
Error Analysis
K
ni;j

m fL
· V
ni;j  1  2rΔtK ni;j 
rΔtL
m
0
Theoretical Analysis
 Fi;j tn  1  rΔtg ml
∂ V

Using the Taylor series expansion, the error of ∂x m ∂zl
X
2


m fL
· W
ni;j  Fi;j tn  1  rΔtg
rΔtL 2 ≤ m  l ≤ 3 is OΔx4  Δz4  caused by the interpola-
m
0 tion formulae presented in Appendix A. Because of the usage


n  rΔtL
2 · W
n  rΔt2 L
2 · W

W

n of the third-order implicit Runge–Kutta method and the trun-
i;j i;j i;j
cated differentiator series method for solving the ODE (8),
X
2
the temporal error, caused by the discretization of the tem-

m · Fi;j tn  1  rΔtg;
frΔtL (17)
poral partial derivative, is OΔt3 . Therefore, we conclude
m
0
that the error introduced by the IRK-DSM is OΔt3 
n

n are defined by Δx4  Δz4 . In other words, the IRK-DSM developed in this
where V

i;j , W
ni;j , and W i;j
article has a fourth-order accuracy in space and third-order
n accuracy in time.
V

i;j
L
· V
ni;j ; (18)

Numerical Error
W
ni;j
V
ni;j  1  2rΔtK ni;j ; (19)
To further illustrate the accuracy of our present method,
in the following discussion we will compare the numerical


n
L
· W
n ; results of the IRK-DSM against other methods, such as the
W i;j i;j (20) conventional FD and the fourth-order LWC methods, with
the exact solution of the 2D acoustic-wave equation for the
and the second-order operator L
2 can be obtained from the
homogeneous medium. Under this consideration, we choose
definition of L
as follows
the following 2D initial problem:
L
2
DiagL2 ; L2 ; L2 
  ∂2u ∂2u 1 ∂2u
1 1 1 1 1 1  2
2 2; (22a)

Diag D; D; D; D; D; D : (21) ∂x 2
∂z c ∂t
ρ ρ ρ ρ ρ ρ


For practical calculations, the implementation of the ex- 2πf0 2πf0
u0; x; z
cos  x cos θ0  z sin θ0 ; (22b)
plicit IRK-DSM is divided into three major steps. The com- c c
putational steps are described next.
1. Computing K ni;j by equation (16).

∂u0; x; z 2πf0 2πf0
(a) Use formulae (A3) to (A9) to compute V

i;j via
n
2πf0 sin  x cos θ0  z sin θ0 ;
∂t c c
equation (18);
n (22c)
(b) Use the obtained result V

i;j in the first step (a) and
apply the similar computational formulae as (A3) where c is the velocity of the plane wave; θ0 , the incident
n
to (A9) to compute L
2 · V

i;j and L
2 · V
ni;j via angle at time t
0; f0 , the frequency. The exact solution
equation (21); of this initial problem is:
(c) Substitute these results obtained in steps (a) and (b) 

into equation (16) to obtain K ni;j . x z
ut; x; z
cos 2πf0 t  cos θ0  sin θ0 : (23)
2. Computing K
ni;j by equation (17). In this step, we first α α
compute W
ni;j using equation (19), and then use the simi-
lar steps as computing K ni;j in the first step to compute
n In this numerical experiment, we choose the number of
K
ni;j . The difference is that these vectors V
ni;j and V

i;j pre-
grid points N
201; the frequency, f0
15 Hz; the veloc-
sented in equation (16) are replaced respectively by the


n when computing K
n . ity, c
4000 m=sec; and θ0
π=4. The relative error for
vectors W
ni;j and W i;j i;j the 2D case is defined by
3. Substitute these obtained results K ni;j and K
ni;j into equa-
n1 
tion (9) to obtain the values of V
i;j at the (n  1)th 1
Er %
PN PN
time level.
i
1 j
1 ut n ; xi ; zj 
2

Note that the fourth terms at the right-hand side of equa- X

N X N 1
2
n
tions (16) and (17) can be easily computed because of the × ui;j  utn ; xi ; zj  2
× 100; (24)
known analytical source function F. i
1 j
1
3344 D. Yang, N. Wang, S. Chen, and G. Song

where uni;j is the numerical solution and utn ; xi ; zj  is the

analytical solution.
Figures 1, 2, and 3 plot the relative errors Er versus time
for different spatial and temporal increments, where three
lines of Er corresponding to the IRK-DSM, the fourth-order
LWC method, and the second-order FDM are shown in a
semilog scale. In these figures, the maximum relative errors
for different cases are listed in Table 1. From these error
curves and Table 1 (Δx
Δz
h), we find that Er in-
creases corresponding to the increase in the temporal and/or
spatial increments for all the three methods. Figures 1, 2, and
3 illustrate that the IRK-DSM achieves the highest numerical
accuracy among all three methods.

Stability Criteria
It is well-known that the temporal increment Δt must be
less than or equal to the Courant limit to keep the numerical
calculation stable. In this section, following the ideas by Figure 2. The relative errors of the IRK-DSM, the fourth-order
Richtmyer and Morton (1967) and Guan and Lu (2006), LWC method, and the second-order FDM measured by Er (formula
we derive the stability criteria of the IRK-DSM for 1D and [24]) are shown in a semilog scale for the 2D initial problem (22).
2D cases. Through a series of mathematical operations, The spatial and the temporal increments are 30 m and 4 × 104 sec,
we obtain the following stability condition for the 1D homo- respectively.
geneous case (see Appendix B)
where c denotes the wave velocity, and Δt and Δx, the time
Δt and space increments, respectively. The maximal value αmax
c ≤ αmax ≈ 0:968 (25)
Δx of the Courant number defined by α
c ΔxΔt
(Dablain, 1986;
or Sei and Symes, 1994) can be found in Appendix B.
For comparison, we also give the stability conditions of
Δx the fourth-order and eighth-order LWC methods (Dablain,
Δt ≤ 0:968 : (26) 1986) in Table 2 for the 1D case. From Table 2 we can see
c

Figure 1. The relative errors of the IRK-DSM, the fourth-order Figure 3. The relative errors of the IRK-DSM, the fourth-order
LWC method, and the second-order FDM measured by Er (formula LWC method, and the second-order FDM measured by Er (formula
[24]) are shown in a semilog scale for the 2D initial problem (22). [24]) are shown in a semilog scale for the 2D initial problem (22).
The spatial and the temporal increments are 20 m and 2 × 104 sec, The spatial and the temporal increments are 40 m and 1 × 103 sec,
respectively. respectively.
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3345

Table 1
Comparison of Maximum Er % for Different Cases and Different Methods
Methods Second-Order FDM Fourth-Order LWC IRK-DSM

Case 1: h
20 m 63.2395 0.899355 0.246253

Δt
2 × 104 sec
Case 2: h
30 m 142.881 4.83214 1.40869
Δt
4 × 104 sec
Case 3: h
40 m 209.775 14.4097 4.55115
Δt
1 × 103 sec

the stability condition of the IRK-DSM is more relaxed than are used (Sei and Symes, 1994; Yang et al., 2002). Therefore,
those of the fourth-order and eighth-order LWC methods. the numerical dispersion is one of the most important issues
For the 2D homogeneous case, the stability condition of which have to be solved in seismic simulations. In this sec-
the IRK-DSM under the condition Δx
Δz
h is given by tion, following the analysis methods presented in the cited
(see Appendix B) references (Dablain, 1986; Yang et al., 2006), we investigate
Δx Δx the numerical dispersion caused by the IRK-DSM through the
Δt ≤ αmax ≈ 0:844 : (27) numerical dispersion relation for the 1D case.
c c
In the following discussion, we show the numerical re-
sults of the IRK-DSM using the dispersion relation (C3) that
These stability criteria for a heterogeneous medium can-
is found in Appendix C. For comparison, we also show nu-
not be directly determined but could be approximated by
merically the dispersion relations of the fourth-order and
using a local homogeneous method. Our conjecture is that
eighth-order LWC methods. The detailed numerical disper-
equations (26) and (27) are approximately correct for a het-
sion analysis of the LWC methods can be found in Yang et al.
erogeneous medium if the maximal value of the wave veloc-
(2006).
ity c is used.
Figures 4, 5, and 6 plot the dispersion relations of the
IRK-DSM (formula C3) and the fourth-order LWC and
Dispersion Analysis and Efficiency eighth-order LWC methods (formula C5 presented in BSSA,
Yang et al., 2006) for the 1D case, which show the ratio R of
Almost all numerical modeling methods suffer from the
the numerical velocity to the real phase velocity versus θ

numerical dispersion when too-coarse grids or too few sam-
kΔx (where k is the wavenumber) for different methods. In
ples per wavelength are used. We may use finer grid to in-
other words, Figures 4–6 show the dispersion errors of the
crease the number of grid points per wavelength or increase
IRK-DSM, the fourth-order LWC, and eighth-order LWC,
the accuracy of numerical method to suppress or eliminate
the numerical dispersion or artifacts. Using fine grid results respectively. From Figure 4 we can see that the numerical
in more computational costs and storages (the total number velocity of the IRK-DSM gradually approximates the exact
of grid points for a fixed-size model will rapidly increase, wave velocity when the Courant number α increases in
especially in 3D cases). It is difficult to apply these tech- the frequency range. It also shows that at the stability limit
Δt
niques in large-scale computation, especially for large-scale (α
c Δx ≈ 0:968), the numerical dispersion of the IRK-
3D simulation of seismic-wave propagation because of its DSM is minimal. This suggests that we can minimize the
intensive use of central processing unit (CPU) time and its numerical dispersion caused by the IRK-DSM by using the
need for large amounts of direct-access memory. Higher- maximum Courant number. Moreover, the change of numer-
order FD or compact FD (e.g., eighth-order method, Dablain, ical dispersion corresponding to the change of the Courant
1986), and staggered-grid FD (Virieux, 1986; Fornberg, number α is not sensitive as compared with the fourth-order
1990) can reduce the numerical dispersion, but merely in- and eighth-order LWC methods. This property is useful in
creasing the accuracy of methods does not result in a propor- practical calculations because we do not need to consider
tional decrease of the numerical dispersion (Yang et al., how to choose an appropriate Courant number in the stability
2006). In other words, the higher-order FD and staggered- range. From Figures 4 and 5 we observe that the numerical
grid FD methods with local operators still suffer from the velocity for the IRK-DSM and the fourth-order LWC method
numerical dispersion when too few samples per wavelength are usually smaller than the exact velocity; this shows that
the numerical dispersion of the IRK-DSM and the fourth-
Table 2 order LWC method follow the exact signal. However, as we
Approximate Maximum Courant Numbers of Different observe in Figure 6, the numerical dispersion of the eighth-
Methods for the 1D Case order LWC method is irregular for a fixed Courant number. It
Methods IRK-DSM OΔt4  Δx4  LWC OΔt4  Δx8  LWC
sometimes leads the exact signal and sometimes follows the
signal. This phenomenon indicates that it is hard to choose
αmax 0.968 0.95 0.929
a suitable Courant number for the eighth-order LWC.
3346 D. Yang, N. Wang, S. Chen, and G. Song

Figure 6. The ratio R of the numerical velocity to the phase

Figure 4. The ratio R of the numerical velocity to the phase
velocity versus wavenumber θ
kΔx for the eighth-order LWC
velocity versus wavenumber θ
kΔx for the IRK-DSM, where
method, where four lines correspond to α
0:3, 0.6, 0.8, and
four lines correspond to α
0:3, 0.6, 0.8, and αmax , respectively.
αmax , respectively.

Comparing Figure 4 with Figures 5 and 6, we can see that for

the same Courant number, the numerical velocity of the IRK- where the source, located at the center of the computation
DSM approximates the exact wave velocity best, especially in domain, is an explosive source and the time variation of
the high-frequency range. It illustrates that the IRK-DMS has the source function is f
sin2πf0  exp4π2 f20 t2 =16
the least numerical dispersion among the three methods. (Zahradnik et al., 1993), and has a Ricker wavelet with
We now investigate the numerical dispersion and com- the peak frequency of f0
25 Hz. The acoustic velocity
putational efficiency of the IRK-DSM through modeling is c0
4000 m= sec, and the computational domain is 0 <
seismic-wave fields for the 2D case. For this case, we con- x ≤ 24 km and 0 < z ≤ 24 km.
sider the following acoustic-wave equation in a 2D homoge- Figure 7 shows the wave-field snapshots at t
2 sec on
neous medium: a coarse grid (Δx
Δz
80 m), generated respectively by
the LWC methods with accuracies of orders 4 and 8, and the

2 
∂ 2u IRK-DSM with the same Courant number α
0:65 and the
2 ∂ u ∂ 2u

c   f; (28) same grid points per wavelength G
f0c·Δx0

2. We can see
∂t2 0
∂x2 ∂z2
that the wave fronts of seismic waves shown in Figure 7,
simulated by the three methods at the same time, are nearly
identical; it took the fourth-order LWC, the eighth-order
LWC, and the IRK-DSM about 6 sec, 8 sec, and 73 sec to
generate Figures 7a,b,c, respectively. This means the compu-
tational cost of the IRK-DSM is more expensive than that of
the high-order LWC methods with the same time steps and
grid points for a fixed-size model. In other words, under the
condition of the same grid interval and time step, our
IRK-DSM method is more expensive than high-order LWC
methods because our method computes more quantities
(not only the wave field itself, but also its time derivative,
the spatial gradients, etc.). However, Figure 7c, generated
by the IRK-DSM, shows that the IRK-DSM has much less nu-
merical dispersion even when the spatial increment is chosen
by Δx
Δz
80 m, whereas the fourth-order LWC and
eighth-order LWC methods exhibit serious numerical disper-
sions (see Figs. 7a,b).
To reduce or eliminate the serious numerical dispersions
Figure 5. The ratio R of the numerical velocity to the phase
velocity versus wavenumber θ
kΔx for the fourth-order LWC generated by the high-order LWC methods, here we set the
method, where four lines correspond to α
0:3, 0.6, 0.8, and Courant number α
0:65 constant and diminish the space
αmax , respectively. increment of the high-order LWC methods, then the grid
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3347

Figure 8. Snapshots of seismic-wave fields at time 2 sec on a

fine grid, generated by (a) the fourth-order LWC method (Δx

Δz
10 m), and (b) the eighth-order LWC method (Δx

Δz
15 m).

number of grid points is 301 × 301 on a coarse grid for

generating Figure 7c. Both of the high-order LWC methods
n1
only need 3 arrays to store the displacement ui;j , uni;j , and
n1
ui;j at each grid point, but the number of grid points on a
Figure 7. Snapshots of seismic-wave fields at time 2 sec on the fine grid for generating Figures 8a and 8b goes up to
coarse grid (Δx
Δz
80 m), generated by (a) the fourth-order
LWC method, (b) the eighth-order LWC method, and (c) the IRK- 2401 × 2401 for the fourth-order LWC and 1601 × 1601
DSM, respectively. for the eighth-order LWC, respectively. It indicates that
the space storage of the IRK-DSM requires only about
12.6% of space storage of the fourth-order LWC method
points per wavelength of the fourth-order LWC and and roughly 28.2% of that of the eighth-order LWC method.
the eighth-order LWC methods increased. Figure 8 shows
the wave-field snapshots at t
2 sec on a fine grid and
under the same Courant number of α
0:65, generated Numerical Modelling
by the fourth-order LWC method (Δx
Δz
10 m, G

c0
f·Δx
20), and the eighth-order LWC method (Δx

In this section we investigate the efficiency of the
Δz
15 m, G
16). IRK-DSM for the 2D isotropic case and the transversely
Comparison between Figure 7c and Figure 8 demon- isotropic model with a vertical symmetry axis (VTI). It
strates that the IRK-DSM with large grid increments should be mentioned in our present numerical experiments
(Δx
Δz
80 m, G
2) can effectively suppress the that wave velocities vary with propagation directions in an-
numerical dispersion, whereas high-order LWC methods isotropic media, so the stability condition Δt ≤ 0:844h=vmax
have to use smaller grid increments to suppress numerical obtained previously is chosen where h (if Δx
Δz
h)
dispersion. The computation time for the IRK-DSM when corresponds to the spatial increment and vmax is the maxi-
mum P-wave or qP-wave velocity in the following models.
larger grids are used is much less than the times by the high-
All our experiments are implemented on a 2-core pentium 4
order LWC methods when smaller grids are used for sup-
computer with 1G memory.
pressing the numerical dispersion. Speaking in detail, in this
experiment example it took the IRK-DSM about 1.22 min
to generate Figure 7c, whereas the high-order LWC methods Isotropic Model
took about 62.77 min and 20.55 min to generate Figures 8a,b,
respectively. It means that the computational speed of the Under this case of our consideration, the elastic-wave
IRK-DSM when the big grid is used is roughly 51 times that equation is
of the fourth-order LWC and about 17 times that of the
∂ 2 u1 ∂ 2 u1 ∂ 2 u1 ∂ 2 u3
eighth-order LWC method on a fine grid to achieve the ρ
λ  2μ  μ  λ  μ  f1
same dispersion accuracy of the IRK-DSM. Meanwhile, ∂t2 ∂x2 ∂z2 ∂x∂z


the storage space required for computation in the IRK-DSM ∂2u ∂ 2 u2 ∂ 2 u2
ρ 22
μ  2  f2
when the big grid is used is also different from that of the ∂t ∂x2 ∂z
high-order LWC methods when smaller grids are used for ∂ 2 u3 ∂ 2 u1 ∂2u ∂2u
suppressing the numerical dispersion. The IRK-DSM needs ρ
λ  μ  μ 23  λ  2μ 23  f3 ;
n1 n1 ∂t 2 ∂x∂z ∂x ∂z
24 arrays to store uni;j , wni;j , ui;j , wi;j , gradients of
(29)
the displacement and the particle velocity, and the middle
variables of K ni;j and L
· V
ni;j at each grid point, and the where λ and μ are the Lamé constants.
3348 D. Yang, N. Wang, S. Chen, and G. Song

In the numerical experiment, we choose the medium The receiver is located at Rx; z
5:85 km; 6:75 km;
constants λ
2:75 GPa and μ
6:25 GPa, and the density the rest of the parameters are the same as those mentioned
ρ
2:1 g=cm3 . The source is an explosive source that is at previously.
the center of the computational domain and has a Ricker Figure 10 presents the waveforms of three components
wavelet with a peak frequency of f0
15 Hz. The time in the isotropic medium at the receiver R, generated by the
variation of the source function is IRK-DSM and the fourth-order LWC, respectively. From
Figure 10 we can see that the IRK-DSM has less or no visible
f1
f2
f3 numerical dispersion, whereas the fourth-order LWC suffers
from the numerical dispersion.

5:76f20 1  160:6f0 t  12  exp80:6f0 t  12 :
(30) VTI Model

To compute the displacement component u2 , we choose Now we consider a two-layered VTI model, for which
space and time steps of Δx
Δz
45 m and Δt
the elastic constants and densities in the upper and lower
0:02 sec, corresponding to the Courant number of layers are listed in Table 3. The spatial and time increments
Δt
α
c Δx
0:767 because of the S-wave velocity of are chosen by Δx
Δz
20 m and Δt
0:0092 sec for
1:725 km=sec. For horizontal component u1 and vertical computing the component u2 , and Δt
0:003 sec for com-
component u3 , we choose spatial and time increments by puting both the horizontal component u1 and the vertical
Δx
Δz
45 m and Δt
0:01 sec, resulting in the max- component u3 , respectively. The computational domain is
imal Courant number of 0.599 because of the P-wave veloc- 0 < x, z ≤ 6 km. The source is located at xs ; zs 

ity of 2:695 km=sec. The computational domain is 0 < x ≤ 3 km; 2:88 km. The rest of the computational parameters
13:5 km and 0 < z ≤ 13:5 km; the number of grid points and the computer environment are the same as those in the
is 301 × 301. isotropic model.
Figure 9 shows the wave-field snapshot at t
1:6 sec, The wave-field snapshots of the three components for
generated by the IRK-DSM. From the three-component snap- u1 , u2 , and u3 at time 0.96 sec are presented in Figure 11.
shots, we can see that the P and SV waves presented in the Figure 11 illustrates that the IRK-DSM has quite less or nearly
horizontal component u1 (Fig. 9a) and vertical component no numerical dispersion and source-generated noises (arti-
u3 (Fig. 9c) are very clear even when a coarse grid facts due to source location at grid points), even for large
(Δx
Δz
45 m) is chosen. It indicates that our present vertical velocity contrast that reaches about 46.4% calcu-
method enables wave propagation to be simulated in large- lated from Table 3. Figure 11 also shows numerous reflected,
scale models through using the coarse computation grids. transmitted, and converted phases. We can see from
For comparison, we present the waveforms computed Figure 11b that the wavefront of the qSH-wave is an ellipse
using the IRK-DSM and the fourth-order LWC method. and the qP- and qSV-waves (Figs. 11a,c) show the direc-
tional dependence on the propagation velocity. The qSV
wavefronts can have cusps and triplications depending on
the value of C13 (Faria and Stoffa, 1994). Triplications
can be observed in the horizontal component qSV-wavefront
(shown in Fig. 11a) and in the vertical component qSV-
wavefront presented in Figure 11c. Moreover, we can
also observe that the arrival times of qSH- and qSV-waves

Figure 9. Snapshots of seismic-wave fields for three compo- Figure 10. A comparison of three-component waveforms in a
nents at time 1.6 sec in an isotropic medium, generated by homogeneous isotropic medium. (a), (b) The synthetic seismograms
the IRK-DSM, for (a) u1 component, (b) u2 component, and (c) u3 are generated by the IRK-DSM and the fourth-order LWC method,
component. respectively.
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3349

Table 3 We also obtain the stability criteria of the IRK-DSM for

Medium Parameters Used in the VTI Model solving 1D and 2D scalar wave equations, and the maximal
Layer
Courant number αmax of the IRK-DSM are 0.968 and 0.844,
(No.) Thickness C11 C13 C33 C44 C66 ρs g cm3 respectively, which are slightly larger than those of the
1 3000 14.2 5.4 18 6.5 3.8 3.2 fourth-order and eighth-order LWC methods as shown in
2 3000 40.8 13.2 50.6 25 13.8 4.2 Table 2 because of combing the implicit RK method with
the truncated differentiator series. Numerical dispersion anal-
ysis shows that the numerical dispersion of the IRK-DSM de-
are different through comparing Figure 11b with Figures 11a creases as the Courant number increases, and the numerical
and c. dispersion reaches its minimum when the Courant number α
equals the maximum Courant number αmax. So we can mini-
mize the numerical dispersion of the IRK-DSM by choosing
Conclusion and Discussion an optimal Courant number that approximates or equals the
In this article, we propose an explicit IRK-DSM via the maximum Courant number αmax in practical seismic-wave
implicit Runge–Kutta method and truncated differentiator modeling. Generally speaking, the numerical dispersion of
series method to solve the acoustic- and elastic-wave equa- the IRK-DSM is smaller than that of the fourth-order and
tions. We first transform the wave equation into a system eighth-order LWC methods (Dablain, 1986) under the same
of ordinary differential equations, and then use an implicit Courant number (see Figs. 4, 5, and 6) especially in the
Runge–Kutta method associated with the differentiator series high-frequency range. So, when the velocity and spatial
method to solve the ODEs. Specifically, we first transform increments are constants, through increasing the temporal in-
the original wave equation (2) into equation (8). Then we crement under keeping the IRK-DSM stable, we can not only
use the interpolation approximations (A3)–(A9) to approxi- lessen the numerical dispersion, but we can also reduce the
mate the high-order spatial derivatives on the right-hand side computational cost and save the storage for computer codes.
of equation (8), converting equation (8) into a system of We have noticed previously that the CPU time of the
semidiscrete ordinary differential equations (ODEs). Finally, IRK-DSM is more expensive than that of the fourth-order
on the basis of both the implicit RK method and the differ- and eighth-order LWC methods (Dablain, 1986) under the
entiator series method, we develop the explicit IRK-DSM condition of the same grid interval and time step. However,
using both the truncated difference-operator Taylor expan- because the IRK-DSM yields less numerical dispersion than
sion and the split-operator method. that of the high-order finite-difference methods such as high-
The error analysis shows that the IRK-DSM is third-order order LWC methods when big grid is used, we can reduce the
accuracy in time and fourth-order accuracy in space; the computational cost and storage space through using a larger
numerical error is less than those of the fourth-order LWC time increment and a coarse spatial grid for a fixed-size
method and second-order FDM under the same condition. model to achieve the same dispersion accuracy as those of
the high-order LWC methods on a finer spatial grid with
smaller time increment for suppressing the numerical disper-
sion. As is confirmed in the numerical dispersion section, the
computational speed of the IRK-DSM when big grid is used is
roughly 51 times that of the fourth-order LWC and about
17 times that of the eighth-order LWC method on a fine grid
to achieve the same dispersion accuracy of the IRK-DSM; the
space storage of the IRK-DSM requires only roughly 12.6%
of the fourth-order LWC method and about 28.2% of that of
the eighth-order LWC method, respectively.
In this IRK-DSM, when determining these high-order
spatial derivatives included in equation (8), the IRK-DSM
uses not only the values of the displacement U and the
particle velocity W at the mesh point (i, j) and its neighbor-
ing grid points (see equations [A3]–[A9]), but also the values
of the gradients of the displacement U and the particle ve-
locity W. Based on such a structure, the IRK-DSM retains
more wave-field information included in the displacement
function, the particle velocity, and their gradients. As a
result, the IRK-DSM can effectively suppress the numerical
Figure 11. Snapshots of seismic-wave fields for three compo- dispersion and source-generated noises caused by discretiz-
nents at time 0.96 sec on the spatial step Δx
Δz
20 m for ing the wave equations when too-coarse grids are used or
(a) u1 component, (b) u2 component, and (c) u3 component. when the models have a large velocity contrast between
3350 D. Yang, N. Wang, S. Chen, and G. Song

adjacent layers. Meanwhile, the IRK-DSM has higher spatial Igel, H., P. Mora, and B. Riollet (1995). Anisotropic wave propagation
accuracy though the method only uses a local difference- through finite-difference grids, Geophysics 60, 1203–1216.
Johnson, C. (1990). Adaptive finite element methods for diffusion
operator to approximate the high-order partial differential op- and convection problems, Comput. Meth. Appl. Mech. Eng. 82,
erator included in wave equations, in which three grid points 301–322.
are used in a spatial direction. Besides, this IRK-DSM pro- Ke, H. L. (1996). Calculation of generalized integrals using differen-
vides more extra seismic information because it simulta- tiator series method, J. Chongqing Jianzhu Univ. 18, 110–113
(in Chinese).
neously computes displacement field, its time derivatives,
Ke, H. L., and H. X. Xie (1999). Differentiator series solution of linear
and their spatial gradients. We initiate possible, more com- differential ordinary equation, Appl. Math. Mech. 20, 59–66.
plicated applications of the IRK-DSM in large-scale seismic Kelly, K., R. Ward, S. Treitel, and R. Alford (1976). Synthetic seismograms:
modeling, reverse time migration, and inversion based on the A finite-difference approach, Geophysics 41, 2–27.
wave equations, for which the computation time and memory Komatitsch, D., and J.-P. Vilotte (1998). The Spectral Element method: An
requirements are the bottleneck for their vast applications. efficient tool to simulate the seismic response of 2D and 3D geological
structures, Bull. Seismol. Soc. Am. 88, 368–392.
Komatitsch, D., C. Barnes, and J. Tromp (2000). Simulation of anisotropic
Data and Resources wave propagation based upon a spectral-element method, Geophysics
65, 1251–1260.
All data used in this article were produced synthetically. Kondoh, Y., Y. Hosaka, and K. Ishii (1994). Kernel optimum nearly ana-
lytical discretization algorithm applied to parabolic and hyperbolic
Acknowledgments equations, Comput. Math. Appl. 27, 59–90.
Kosloff, D., and E. Baysal (1982). Forward modeling by a Fourier method,
We thank Fred F. Pollitz (associate editor) and one anonymous Geophysics 47, 1402–1412.
reviewer for their detailed comments and helpful suggestions that greatly Kosloff, D., M. Reshef, and D. Loewenthal (1984). Elastic wave calculations
contributed to improving the article. Special thanks to Biaolong Hua and by the Fourier method, Bull. Seismol. Soc. Am. 74, 875–891.
Jiming Peng for their revisions that substantially improved the presentation Lele, S. K. (1992). Compact finite difference scheme with spectral-like
of the article. This work was supported by the National Science Fund for resolution, J. Comp. Phys. 103,16–42.
Distinguished Young Scholars of China (Grant No. 40725012) and the Levander, A. R. (1988). Fourth-order finite-difference P-SV seismograms,
National Major Fundamental Research and Development Project of China Geophysics 53, 1425–1436.
(No. 2007CB411704). Priolo, E., and G. Seriani (1991). A numerical investigation of Chebyshev
spectral element method for acoustic wave propagation, Proc. of the
References 13th World Congress on Computation and Applied Mathematics,
Vol. 2, 551–556.
Alford, R. M., K. R. Kelly, and D. M. Boore (1974). Accuracy of finite- Priolo, E., J. M. Carcione, and G. Seriani (1994). Numerical simulation of
difference modeling of the acoustic wave equation, Geophysics 39, interface waves by high-order spectral modeling techniques, J. Acoust.
834–842. Soc. Am. 95, 681–693.
Blanch, J. O., and J. O. A. Robertsson (1997). A modified Lax–Wendroff Richtmyer, R. D., and K. W. Morton (1967). Difference Methods for Initial
correction for wave propagation in media described by Zener elements, Value Problems, Interscience, New York.
Geophys. J. Int. 131, 381–386. Sei, A., and W. Symes (1994). Dispersion analysis of numerical wave
Booth, D. C., and S. Crampin (1983a). The anisotropic reflectivity tech- propagation and its computational consequences, J. Sci. Comput.
nique: Theory, Geophys. J. R. Astr. Soc. 72, 755–765. 10, 1–27.
Booth, D. C., and S. Crampin (1983b). The anisotropic reflectivity tech- Takeuchi, N., and R. J. Geller (2000). Optimally accurate second order time-
nique: Anomalous arrives from an anisotropic upper mantle, Geophys. domain finite difference scheme for computing synthetic seismograms
J. R. Astr. Soc. 72, 767–782. in 2D and 3D media, Phys. Earth Planet. Interiors 119, 99–131.
Dablain, M. A. (1986). The application of high-order differencing to the Turner, M. J., R. W. Clough, H. C. Martin, and L. J. Topp (1956). Stiffness
scalar wave equation, Geophysics 51, 54–66. and deflection analysis of complex structures, J. Aeronaut. Sci. 23,
Eriksson, K., and C. Johnson (1991). Adaptive finite element methods for 805–823.
parabolic problems I: A linear model problem, SIAM J. Numer. Anal. Virieux, J. (1986). P-SV wave propagation in heterogeneous media:
28, 43–77. Velocity-stress finite-difference method, Geophysics 51, 889–901.
Faria, E. L., and P. L. Stoffa (1994). Finite-difference modelling in trans- Vlastos, S., E. Liu, I. G. Main, and X. Y. Li (2003). Numerical simulation
versely isotropic media, Geophysics 59, 282–289. of wave propagation in media with discrete distributions of fractures:
Fei, T., and K. Larner (1995). Elimination of numerical dispersion in finite- Effects of fracture sizes and spatial distributions, Geophys. J. Int. 152,
difference modeling and migration by flux-corrected transport, Geo- 649–668.
physics 60, 1830–1842. Wang, S. Q., D. H. Yang, and K. D. Yang (2002). Compact finite difference
Fornberg, B. (1990). High-order finite differences and pseudo-spectral scheme for elastic equations, J. Tsinghua Univ. (Sci. & Tech.) 42,
method on staggered grids, SIAM J. Numer. Anal. 27, 904–918. 1128–1131 (in Chinese).
Fuchs, K., and G. Müller (1971). Computation of synthetic seismograms Whiteman, J. R. (1975). A Bibliography for Finite Elements, Academic
with the reflectivity method and comparison with observations, Press, New York.
Geophys. J. Int. 23, 417–433. Yang, D. H., E. Liu, and Z. J. Zhang (2008). Evaluation of the u  W finite
Geller, R. J., and N. Takeuchi (1998). Optimally accurate second-order time- element method in anisotropic porous media, J. Seism. Expl. 17,
domain finite difference scheme for the elastic equation of motion: 273–299.
One-dimensional case, Geophys. J. Int. 135, 48–62. Yang, D. H., E. Liu, Z. J. Zhang, and J. Teng (2002). Finite-difference mod-
Guan, Z., and J. F. Lu (2006). Numerical Methods, Tsinghua University eling in two-dimensional anisotropic media using a flux-corrected
Press, Beijing (in Chinese). transport technique, Geophys. J. Int. 148, 320–328.
Hairer, E., S. P. Nørsett, and G. Wanner (1993). Solving Ordinary Yang, D. H., J. M. Peng, M. Lu, and T. Terlaky (2006). Optimal nearly
Differential Equations I: Nonstiff Problems, Springer-Verlag, Berlin, analytic discrete approximation to the scalar wave equation, Bull.
Heidelberg. Seismol. Soc. Am. 96, 1114–1130.
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3351

1 1
Yang, D. H., G. J. Song, S. Chen, and B. Y. Hou (2007). An improved nearly
∂xz Uni;j
E1  E1 n
x ∂z U i;j  E1  E1 n
z ∂x Ui;j
analytical discrete method: An efficient tool to simulate the seismic 2Δx x 2Δz z
response of 2D porous structures, J. Geophys. Eng. 4, 40–52.
1
Yang, D. H., J. W. Teng, Z. J. Zhang, and E. Liu (2003). A nearly-analytic  E1 E1  E1x E1 1 1 1 1 n
z  Ex Ez  Ex Ez Ui;j ;
discrete method for acoustic and elastic wave equations in anisotropic 4ΔxΔz x z
media, Bull. Seismol. Soc. Am. 93, 882–890. (A5)
Zahradnik, J., P. Moczo, and T. Hron (1993). Testing four elastic finite-
difference schemes for behavior at discontinuities, Bull. Seismol.
15
Soc. Am. 83, 107–129. ∂3x V ni;j
E1x  E1 n
x V i;j
Zhang, Z. J., G. J. Wang, and J. M. Harris (1999). Multi-component wave- 2Δx3
field simulation in viscous extensively dilatancy anisotropic media, 3
Phys. Earth Planet. Interiors 114, 25–38.  E1x  8I  E1 n
x ∂x V i;j ; (A6)
2Δx2
Zheng, H. S., Z. J. Zhang, and E. R. Liu (2006). Nonlinear seismic wave
propagation in anisotropic media using the flux-corrected transport
15
technique, Geophys. J. Int. 165, 943–956. ∂3z V ni;j
E1z  E1 n
z V i;j
2Δz3
3
Appendix A  E1z  8I  E1 n
z ∂z V i;j ; (A7)
2Δz2
1
∂2xz V ni;j
5E1x E1z  5E1 1 1 1 1 1
x Ez  Ex Ez  Ex Ez
Evaluation of High-Order Derivatives 4Δx2 Δz
 4E1z  4E1 1 n
z  6Ex  6Ex V i;j
1
When we use the explicit IRK-DSM to compute the
values of U at time tn1 in synthetic seismograms, we first 1
 E1x E1z  E1 1
x Ez  Ex  Ex
1 1
need the high-order space derivatives included in equa- 2ΔxΔz
tions (9), (16), and (17), or equation (8). To determine the 1 2
 2δ2z ∂x V ni;j  δ x ∂z V ni;j ; (A8)
high-order derivatives, following Konddoh et al. (1994) Δx2
and Yang et al. (2003), we introduce the interpolation func-
tion of the spatial steps Δx and Δz as follows 1
∂x2z V ni;j
5E1x E1z  5E1 1 1 1 1 1
x Ez  Ex Ez  Ex Ez

 4Δx2 Δz
X5
1 ∂ ∂ r
GΔx; Δz
Δx  Δz V; (A1)  4E1x  4E1 1 1 n
x  6Ez  6Ez V i;j
r
0
r! ∂x ∂z
1
 E1x E1z  E1 1
x Ez  Ez  Ez
1 1
and construct the connection relations between the point 2ΔxΔz
(i, j) and its neighboring nodes. For example, at the grid 1
point (i  1, j) we have the following connection relations,  2δ2x ∂z V ni;j  2 δ 2z ∂x V ni;j ; (A9)
Δz
GΔx; 0ni;j
V ni1;j ; where δ2z V ni;j
V ni;j1  2V ni;j  V ni;j1 , E1z V ni;j
V ni;j1 ,
 n
n
∂ ∂ E1 n n
z V i;j
V i;j1 . The operators δ x , Ex , and Ex
2 1 1
can be
GΔx; 0
V ; similarly defined. And V i;j , ∂x V i;j , ∂z V i;j , and ∂mxkz V ni;j
n n n
∂Δx i;j ∂x i1;j ∂ ∂
 n
n denote ViΔx; jΔz; nΔt, ∂x ViΔx; jΔz; nΔt, ∂z ViΔx;
∂ ∂ jΔz; nΔt, and ∂ mk m k n
V=∂x ∂z i;j , respectively.
GΔx; 0
V (A2)
∂Δz i;j ∂z i1;j

and similarly the remaining 21 connection relations at the Appendix B

other 7 neighboring nodes can be easily written. Notice that
the vector V is defined by V
U; WT, where U and W are Derivation of Stability Criteria
the displacement and the particle velocity, respectively.
From the 24 relations at the nodes (i  1, j), (i  1, j),
(i, j  1), (i, j  1), (i  1, j  1), (i  1, j  1), (i  1, 1D Homogeneous Case
j  1), and (i  1, j  1), we can obtain the following ap- For the 1D case, these formulae (A3) and (A6) can be
proximation formulae, which are the similar expressions to degenerated into the following formulae:
those suggested by Yang et al. (2003, 2007). For conve-
nience, we present the approximation formulae as follows 2 2 n 1
∂2x V nj
δx V j  E1  E1 n
x ∂x V j ; (B1)
Δx2 2Δx x
2 2 n 1
∂2x V ni;j
δ x V i;j  E1  E1 n
x ∂x V i;j ; (A3)
Δx2 2Δx x
15
∂3x V nj
E1x  E1 n
x V j
2Δx3
2 2 n 1 3
∂2z V ni;j
δ V 
2 z i;j
E1  E1 n
z ∂z V i;j ; (A4)  E1x  8I  E1 n
x ∂x V j : (B2)
Δz 2Δz z 2Δx2
3352 D. Yang, N. Wang, S. Chen, and G. Song

To obtain the stability condition of the IRK-DSM, we consid- g21
e3iθ α2 fα4 r3 1:375  2:75r
er the harmonic solution of equations (9), (16), and (17) for
the 1D case. Substituting the following solution  α4 e6iθ 1:375  2:75rr3
0 n 1  α2 eiθ r0:25  0:25r  27α2 r2  54α2 r3 
u
B w n C  α2 e5iθ r0:25  0:25r  27α2 r2  54α2 r3 
V
nj
B C
@ ∂x un A expikjh; (B3)
 e2iθ 2  α4 129:375  258:75rr3
n
∂x w
 α2 r16  16r  e4iθ 2  α4 129:375
into equations (9), (16), and (17) with relations (B1) and  258:75rr3  α2 r16  16r
(B2), we can obtain the following equation
 e3iθ 4  α2 31:5
V
n1
GV
n ; (B4)  31:5rr  α4 r3 310  620rg=Δt (B10)

where the amplification matrix G is

2 3 g22
1  α2 2  eiθ  eiθ   α4 e2iθ r2 0:375
g11 g12 g13 g14
6 g21 g22 g23 g24 7  e2iθ 47:25  94:5r  e4iθ 0:375  0:75r
G
6 4 g31 g32 g33
7 (B5)
g34 5
 0:75r  eiθ 24  48r  e3iθ 24  48r
g41 g42 g43 g44
 α6 e3iθ r4 0:6875  e2iθ 64:6875  129:375r
with
 e4iθ 64:6875  129:375r  eiθ 13:5  27r
iθ iθ 4 2iθ 2
g11
1  α 2  e2
e α e r 0:375
 e5iθ 13:5  27r  1:375r
 e 47:25  94:5r  e 0:375  0:75r  0:75r
2iθ 4iθ
 e6iθ 0:6875  1:375r
 eiθ 24  48r  e3iθ 24  48r
 e3iθ 155  310r (B11)
 α6 e3iθ r4 0:6875  e2iθ 64:6875  129:375r
 e4iθ 64:6875  129:375r  eiθ 13:5  27r

 e5iθ 13:5  27r  1:375r 1
g23
α4 h 8eiθ  8eiθ  0:25e2iθ  0:25e2iθ r
Δt
 e6iθ 0:6875  1:375r  e3iθ 155  310r 
1 iθ iθ 2iθ
(B6)  8e  8e  0:25e  0:25e r 
2iθ 2
Δt

g12
e2iθ fα4 0:25  0:5rr3  α4 e4iθ 0:25  0:5rr3  α6 h
1
111:25eiθ  111:258eiθ  e2iθ  e2iθ
Δt
 α2 eiθ r2  2r  16α2 r2  32α2 r3 
1
 α2 e3iθ r2  2r  16α2 r2  32α2 r3   0:25e3iθ  0:25e3iθ r3   222:5eiθ
Δt

 e2iθ 1  α4 31:5  63rr3  α2 r4  4rgΔt iθ 2iθ 3iθ
 222:5e  2e  2e 2iθ
 0:5e  0:5e r 
3iθ 4
(B7)
α2 eiθ 0:5  0:5e2iθ h
iθ
g13
α 0:25e h  0:25e
2 iθ
h  α h12e
4 iθ  (B12)
Δt
 12eiθ  0:375e2iθ  0:375e2iθ r2  24eiθ
 24eiθ  0:75e2iθ  0:75e2iθ r3   α6 h55:625eiθ
g24
α2 0:25eiθ h  0:25eiθ h
 55:625eiθ  0:5e2iθ  0:5e2iθ  0:125e3iθ
 α4 h12eiθ  12eiθ
 0:125e3iθ r4  111:25eiθ  111:25eiθ
 0:375e2iθ  0:375e2iθ r2
 e2iθ  e2iθ  0:25e3iθ  0:25e3iθ r5  (B8)
 24eiθ  24eiθ  0:75e2iθ  0:75e2iθ r3 
 α6 h55:625eiθ  55:625eiθ  0:5e2iθ
g14
α2 eiθ h0:5  e2iθ 0:5  0:5rrΔt  α4 h8eiθ  0:5e2iθ  0:125e3iθ  0:125e3iθ r4
 8eiθ  0:25e2iθ  0:25e2iθ r3 Δt  111:25eiθ  111:25eiθ  e2iθ
 16eiθ  16eiθ  0:5e2iθ  0:5e2iθ r4 Δt (B9)  e2iθ  0:25e3iθ  0:25e3iθ r5  (B13)
An Explicit Method Based on the Implicit Runge–Kutta Algorithm for Solving Wave Equations 3353



3:75eiθ  3:75eiθ 1 3:75eiθ  3:75eiθ 1
g31
α 2
 fα4 180eiθ g42
α2  fα4 180eiθ
h h h h
 180eiθ  5:625e2iθ  5:625e2iθ r2  180eiθ  5:625e2iθ  5:625e2iθ r2

 360eiθ  360eiθ  11:25e2iθ  11:25e2iθ r3 g  360eiθ  360eiθ  11:25e2iθ  11:25e2iθ r3 g
1 1
 fα6 834:375eiθ  834:375eiθ  7:5e2iθ  fα6 834:375eiθ  834:375eiθ  7:5e2iθ
h h
 7:5e2iθ  1:875e3iθ  1:875e3iθ r4  7:5e2iθ  1:875e3iθ  1:875e3iθ r4

 1668:75eiθ  1668:75eiθ  15e2iθ  15e2iθ  1668:75eiθ  1668:75eiθ  15e2iθ  15e2iθ

 3:75e3iθ  3:75e3iθ r5 g (B14)  3:75e3iθ  3:75e3iθ r5 g (B19)

eiθ r7:5  e2iθ 7:5  7:5r  7:5rΔt 1 3iθ 2 4 3

g32
α2 g43
e α fα r 0:375  0:75r
h Δt
1 4
 fα 120eiθ  120eiθ  3:75e2iθ  α4 e6iθ 0:375  0:75rr3
h
 3:75e2iθ r3 Δt  240eiθ  240eiθ  7:5e2iθ  α2 eiθ r1:5  1:5r  24α2 r2  48α2 r3 

 7:5e2iθ r4 Δtg (B15)  α2 e5iθ r1:5  1:5r  24α2 r2  48α2 r3 

 e2iθ 1:5  α2 36  36rr
 α4 r3 661:875  1323:75r
g33
1  α2 6  0:75eiθ  0:75eiθ 
 e4iθ 1:5  α2 36  36rr
 α4 e2iθ r2 2:25  e2iθ 234  468r
 α4 r3 661:875  1323:75r
 eiθ 54  108r  4:5r  e3iθ 54  108r
 e3iθ 12  α2 156  156rr
 e4iθ 2:25  4:5r  α6 e3iθ r4 0:1875
 α4 r3 2100  4200rg (B20)
 eiθ 12  24r  e5iθ 12  24r
 e6iθ 0:1875  0:375r  0:375r
 e2iθ 330:938  661:875r g44
1  α2 6  0:75eiθ  0:75eiθ 
 e3iθ 1050  2100r  α4 e2iθ r2 2:25  e2iθ 234  468r
 e4iθ 330:938  661:875r (B16)  eiθ 54  108r  4:5r  e3iθ 54  108r
 e4iθ 2:25  4:5r  α6 e3iθ r4 0:1875
g34
e2iθ fα4 r3 1:5  3r  α4 e4iθ r3 1:5  3r  eiθ 12  24r  e5iθ 12  24r
 α2 eiθ r1:5  1:5r  36α2 r2  72α2 r3   e6iθ 0:1875  0:375r  0:375r
 α2 e3iθ r1:5  1:5r  36α2 r2  72α2 r3   e2iθ 330:938  661:875r
 e2iθ 1  α4 156  312rr3  e3iθ 1050  2100r
 α2 r12  12rgΔt (B17)  e4iθ 330:938  661:875r (B21)
p
in which θ
kh, h
Δx, r
12  3
6 ,
the Courant number
α4 (Dablain, 1986; Sei and Symes, 1994) is defined by
g41
120eiθ  120eiθ  3:75e2iθ  3:75e2iθ r α
c0 Δt=h with h
Δx.
hΔt
Let G denote the conjugate transpose matrix of G. Fol-
 120eiθ  120eiθ  3:75e2iθ  3:75e2iθ r2 
lowing the analyses of Richtmyer and Morton (1967) and
α6 Guan and Lu (2006), we know that the IRK-DSM with the
 1668:75eiθ  1668:75eiθ  15e2iθ
hΔt amplification matrix G is stable if the spectral radius ρG ·
 15e2iθ  3:75e3iθ  3:75e3iθ r3  3337:5eiθ G of the matrix G · G satisfies ρG · G ≤ 1. So we can
obtain the stability condition, deriving from the condition of
 3337:5eiθ  30e2iθ  30e2iθ  7:5e3iθ ρG · G ≤ 1, as follows
α2 eiθ 7:5  7:5e2iθ 
 7:5e3iθ r4   (B18) α ≤ αmax ≤ 0:968; (B22)
hΔt
 3354 D. Yang, N. Wang, S. Chen, and G. Song

0 1
or u0
B w0 C
h h V
nj
B C
@ ∂x u0 A expiωnum nΔt  kjh (C1)
Δt ≤ αmax ≈ 0:968 ; (B23)
c0 c0 ∂x w0

where αmax denotes the maximum Courant number that into equations (9), (16), and (17) with relations (B1) and (B2)
keeps the numerical calculation stable. to obtain the dispersion equation as follows:

DetM
0; (C2)

where k is the wavenumber. Owing to the complexity of

2D Homogeneous Case elements of the matrix M, we omit the detail expressions
For the 2D problem, we consider the case of Δx
of the matrix M.
Δz
h. Following the same steps as discussed in the 1D From the dispersion relation (C2), we can obtain the
case, we can obtain the following stability condition ratio of the numerical velocity to the phase velocity c0 as
follows:
h h
Δt ≤ αmax ≈ 0:844 : (B24) cnum ωnum Δt γ
c0 c0 R


; (C3)
c0 αθ αθ

where α is the Courant number, θ
kh, h
Δx, and γ

Appendix C ωnum Δt satisfy the dispersion equation (C3).

Derivation of the Numerical Dispersion Relation Department of Mathematical Sciences

Tsinghua University
To minimize the dispersion error, we derive the numer- Beijing 100084, China
ical dispersion relation of the IRK-DSM for the 1D case. For [email protected]
this, following the dispersion analysis methods presented in
Dablain (1986) and Yang et al. (2006), we consider the har-
monic solution of equation (9) and substitute the solution Manuscript received 30 November 2008

View publication stats