IMPLICIT SOLUTIONS TO THE WAVE EQUATION BASED ON THE METHOD OF
LINES TRANSPOSE
By
Gerard Lee Van Groningen

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mathematics
2012

ABSTRACT
IMPLICIT SOLUTIONS TO THE WAVE EQUATION BASED ON THE
METHOD OF LINES TRANSPOSE
By
Gerard Lee Van Groningen
We present a numerical method for computing the wave equation implicitly. The approach discretizes the wave equation in time using the method of lines transpose, also known
as the transverse method of lines or Rothe’s method. This differs from conventional methods
in that we solve the resulting system of ODEs using boundary integral methods. We then
analyze the fully discretized solution resulting from the midpoint and trapezoidal quadrature
rules and show that convergent and unconditionally stable schemes result. We also show that
the choice of discretization in time can lead to various schemes of prescribed accuracy, which
may or may not introduce numerical dissipation into the approximate solution.
We start with the simplest case of solving the wave equation in one dimension using
either a free space solution or Dirichlet or Neumann Boundary conditions. We then analyze
the stability and consistency of the method, as well as investigating the dispersion relations
and deriving the phase error. Next, some numerical examples are presented which give
validation to the error estimates. Further, the method is adapted for both outflow boundary
conditions, using either one-way waves or a perfectly matched absorbing layer, as well as the
implementation of a soft source. Finally, we utilize an ADI scheme to explore solving the
wave equation in higher dimensions.
Since this method of lines transpose approach is implicit, it removes the usual CFL
stability limit inherent in explicit time stepping methods for solving the wave equation, and
thus the algorithm will be more efficient than these explicit methods.

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

v

. . . . . . . . . . . . . . . . . .

vi

. . . . .
. . . . .
. . . . .
Toeplitz
. . . . .

. . . . . .
. . . . . .
. . . . . .
Matrices
. . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

1
1
8
8
12

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

15
15
15
16
17
22

3 Analysis of the Free Space Solution . . . . . .
3.1 Convergence . . . . . . . . . . . . . . . . . . .
3.1.1 Consistency of the MOLT method . . .
3.1.1.1 Non-dissipative Scheme . . .
3.1.1.2 Dissipative Schemes . . . . .
3.1.1.3 Fully Discrete Schemes . . . .
3.1.2 Stability of the MOLT method . . . .
3.1.2.1 The Eigenvalues of Ap . . . .
3.1.2.2 Dissipative Schemes . . . . .
3.1.2.3 Non-dissipative Scheme . . .
3.2 Dispersion . . . . . . . . . . . . . . . . . . . .
3.2.1 Semi-discrete Non-dissipative Scheme .
3.2.2 Semi-discrete Dissipative Schemes . . .
3.2.3 Fully Discrete Non-dissipative Scheme
3.2.4 Fully Discrete Dissipative Schemes . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

27
27
28
28
30
33
36
38
42
45
46
47
50
54
59

List of Figures

. . . . . . . . . . . . . .

1 Introduction . . . . . . . . . . . .
1.1 Literature Review . . . . . . . .
1.2 Prerequisites . . . . . . . . . . .
1.2.1 Eigenvalues of Perturbed
1.2.2 Schur Polynomial Test .

. . . . . . .
. . . . . . .
. . . . . . .
Tri-diagonal
. . . . . . .

2 MOLT Formulation and Integral Solution
2.1 Formulation . . . . . . . . . . . . . . . . .
2.1.1 Dissipative Schemes . . . . . . . . .
2.1.2 Purely Dispersive Schemes . . . . .
2.2 Integral Solution . . . . . . . . . . . . . .
2.3 Numerical Approach . . . . . . . . . . . .

iii

.
.
.
.
.
.

4 Analysis of the Solution on a Bounded Domain
4.1 Consistency . . . . . . . . . . . . . . . . . . . . .
4.1.1 Non-dissipative Scheme . . . . . . . . . . .
4.1.2 Dissipative Schemes . . . . . . . . . . . . .
4.2 Stability . . . . . . . . . . . . . . . . . . . . . . .
4.3 Numerical Results . . . . . . . . . . . . . . . . . .
4.4 Fast Implementation . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

64
64
64
71
72
83
88

5 Non-Reflecting Boundary Conditions .
5.1 Perfectly Matched Layers . . . . . . . .
5.1.1 Constant Damping Function . .
5.1.1.1 Implementation . . . .
5.1.1.2 Implicit Update . . . .
5.1.1.3 Performance . . . . .
5.1.2 Polynomial Damping Functions
5.1.2.1 Implementation . . . .
5.1.2.2 Optimization . . . . .
5.2 One-Way Waves . . . . . . . . . . . . .
5.2.1 Performance . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

90
92
95
95
99
106
108
109
129
130
139

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

6 MOLT in Higher Dimensions Utilizing an Alternating Directly Implicit
Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
6.1 Freespace Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.2 Finite Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
7 Soft Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
7.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
7.2 Higher Dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
8.1 Further Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
8.2 Numerical Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
Bibliography . . . . . . . . . . . . . . .

iv

. . . . . . . . . . . . . . . . . . . 159

LIST OF TABLES

Table 5.1

Reflection and Damping for the Dissipative 1st Order Scheme with a
Piecewise Constant Damping Function . . . . . . . . . . . . . . . . . 108

Table 5.2

Reflection and Damping for the Dissipative 2nd Order Scheme with
a Piecewise Constant Damping Function . . . . . . . . . . . . . . . 108

Table 5.3

Optimization of the Dissipative 1st Order Quadratic Damping Function130

Table 5.4

Performance of the Dissipative 1st Order Scheme with a Non-reflecting
Boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

v

LIST OF FIGURES

Figure 3.1

Relative Phase Error for the Non-dissipative Scheme . . . . . . . . .

60

Figure 4.1

Contour Plot of the Minimum Eigenvalue of the Dissipative Schemes
for Zero Dirichlet Boundary Conditions with Midpoint Quadrature .

84

Contour Plot of the Minimum Eigenvalue of the Dissipative Schemes
for Zero Neumann Boundary Conditions with Midpoint Quadrature

84

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Convergence of the First Order Dissipative Scheme with ν = ∆x/(c∆t) =
0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
√
Convergence of the Second Order Dissipative Scheme with ν/ 2 =
∆x/(c∆t) = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

√
Convergence of the Second Order Non-dissipative Scheme with ν/ 2 =
∆x/(c∆t) = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

Figure 5.1

Damping Profiles for Low Order Polynomials . . . . . . . . . . . . . 115

Figure 6.1

Propagation of the First Order Dissipative ADI Scheme at t = 0 . . 146

Figure 6.2

Propagation of the First Order Dissipative ADI Scheme at t = 0.5 . 147

Figure 6.3

Propagation of the First Order Dissipative ADI Scheme at t = 1.25

Figure 7.1

Wave Introduction via Soft Source . . . . . . . . . . . . . . . . . . . 152

Figure 7.2

Long Term Time Behavior of the Soft Source . . . . . . . . . . . . . 153

vi

147

Chapter 1
Introduction

1.1

Literature Review

The wave equation is an important hyperbolic Partial Differential Equation (PDE) that
arises in acoustics, electromagnetics, and fluid dynamics. The focus of this dissertation is
computing numerical solutions to the wave equation

vtt = c2 2 v,

x ∈ Ω,

v(x, 0) = f (x),

t ≥ 0,

(1.1)

vt (x, 0) = g(x),

where Ω may either consist of the entire region of Rn , n ∈ {1, 2, 3}, or of a subset of Rn , in
which case boundary conditions must be applied.
We are particularly interested in the setting of Maxwell’s equations as they arise in nonrelativistic plasma situations. Maxwell’s equations may be expressed in terms of a scalar

1

potential Ψ and a vector potential A as

E =− Ψ−

∂A
,
∂t

B=

× A,

(1.2)

where E represents the electric field and B represents the magnetic field. If one chooses the
Lorenz gauge condition

1 ∂Ψ
·A+ 2
= 0,
c ∂t

(1.2) may be expressed in the form of a wave equation as
1 ∂ 2Ψ
−
c2 ∂t2

2Ψ

1 ∂ 2A
−
c2 ∂t

ρ
= ,

2A

= µJ,

where ρ is the charge density, is the free space permittivity, µ is the free space permeability
and J is the current. In this situation, we wish to develop a numerical solver which will be
able to handle the time scale gap between the electron plasma frequency and the speed of
light.
There currently exist several well established approaches to computing a numerical solution to the wave equation. One such approach is to use a separable solution. That is, we
assume the ansatz

v(x, t) = ψ(x)e−iωt ,

2

which leads to the oscillatory Helmholtz equation

2ψ

+ k 2 ψ = 0,

where k = ω/c is the wave number and ω is the frequency of an eigenmode. While the
resulting Helmholtz equation can be solved efficiently using fast summation methods [22, 41],
this approach is only useful for generating time harmonic solutions to the wave equation.
A fairly standard approach for hyperbolic PDEs is the Method of Lines (MOL) discretization. The prevalent idea behind the Method of Lines was introduced in 1949 by Courant
and Lax [28] as an adaptation for using the method of characteristics as a numerical solver
for hyperbolic PDEs. Courant and Lax began with semi-linear problems in which the characteristics are known a priori and the solution to the PDE is then found using an iterative
method over the resulting integral equations. However, when this approach was extended
to non-linear problems, problems arose since the characteristics were not known beforehand.
Therefore, the resulting integral equations might no longer be easily solved using the same
iterative methods. Courant and Lax overcame this by replacing the characteristics with
lines upon which the spatial variable x is held constant. This gave rise to a PDE which is
discretized only in space, resulting in a system of Ordinary Differential Equations (ODEs),
which can then be solved using whichever Initial Value Problem (IVP) solver is desired. Each
individual ODE then represents a solution at some point on the grid as a function of time.
Many studies have been done on the numerical properties of the MOL discretization.
Courant et al. [27] studied convergence based on the ratio of step sizes in different dimensions. Peter Lax [52] used his MOL discretization design to prove the existence of solutions
to non-analytic initial value problems that result from quasi-linear hyperbolic PDEs. In
3

addition to its original purpose in solving hyperbolic PDEs, the method of lines has also
been studied for both parabolic [71] and elliptic [49] equations as well. Much work has been
done in establishing the accuracy and stability of the MOL discretization, including Verwer
and Sanz-Serna [69], who showed convergence of both the space discretization based on a
logarithmic matrix norm and convergence of the time integration based on C-stability, as
well as Reddy and Trefethen [60], who used -pseudo-eigenvalues to obtain conditions that
were both necessary and sufficient for the stability of the method of lines.
Another technique whose roots trace back to the method of characteristics are the semiLagrangian methods, first introduced by Courant et al. [27]. In a purely Lagrangian solver
for the advection equation, the exact value of the solution at a given time is found by
tracing back along its characteristic to find the value at the previous time step. However,
this requires that the mesh be flexible and be able to move in time. As this is frequently
impractical to implement, semi-Lagrangian methods require a fixed mesh and then trace back
the characteristic to find its value at the previous time step. In the case of more complicated
PDEs where the characteristic may not be easily found, the characteristic must instead be
numerically approximated. In the likely event that this characteristic does not lie on a mesh
point at the previous time step, polynomial interpolation will be used to determine the value
of the point at which the characteristic does cross. This essentially reduces the problem
to solving a set of independent ODEs which will approximate the characteristic and then
determining the function value at those points using interpolation.
While originally based on advection equations (e.g. [62]), semi-Lagrangian methods are
now very commonly used for a range of problems, such as the Vlasov equation [65], but
especially in the atmospheric sciences [63, 66] due to its beneficial numerical properties, such

4

as the lack of restriction based on a Courant number [19]. Since the actual semi-Lagrangian
method can vary in order based on both the approximation of the characteristic and the order
of the interpolation used, a variety of studies have been done on the stability, convergence
and cost effectiveness of these various implementations [10, 19, 34, 53, 66].
Another common solver for hyperbolic PDEs is the finite volume methods, beginning in
1959 with Godunov’s method [38]. Godunov divided his domain into cells and found the
average of the numerical solution over each cell. He then constructed a Riemann problem
which consisted of the PDE whose initial conditions were given by the piecewise constant
function corresponding to the average numerical value over each cell. Next, he evolved the
solution to this Riemann problem for a time step and then recalculated the average value over
each cell. This process was then repeated until the desired final time was reached. The most
significant benefit of Godunov type schemes was a greater ability to handle shocks compared
to other contemporary numerical solvers, due to procuring the average value of the numerical
solution over each cell. The further study of the convergence and other numerical properties
of Godunov-type schemes have been studied for a range of problems [29, 45, 54, 59]. Current
examples of finite volume methods whose roots trace back to Godunov-type schemes include
the weighted essentially non-oscillatory (WENO) method, which uses an adaptive stencil
and a combination of lower order reconstructions to obtain a higher order solution to the
PDE, as well as the discontinuous Galerkin (DG) method, which constructs a solution over
each element using lower order polynomials before updating the inter-element terms.
Godunov-type schemes are not without their drawbacks, however. Due to the Riemann
problem’s initial condition of a piecewise constant function, the scheme was limited to first
order. This order may be made higher though, if the initial condition uses piecewise poly-

5

nomials instead of piecewise constant functions. However, making this change will increase
the complexity of the problem and higher order schemes may still be hard to achieve. Additionally, the scheme is restricted in that the size of the time step is bounded by half the size
of the spatial step multiplied by the fastest signal in the solution.
Integral equation methods, such as [6, 39, 40, 51], may also be used to solve hyperbolic
problems. The benefits of integral equation methods, in comparison to traditional finite
difference and finite element methods, include a much simpler incorporation of complex
boundary geometries and the easier attainment of higher order. The biggest drawback to
integral equation methods, however, is that they possess a large computational complexity,
with a time independent problem requiring a dense N 2 × N 2 matrix for an N × N grid.
Advances are being made in reducing the amount of necessary operations, most notably the
fast multipole method (FMM) [25, 61] for the wave scattering problem. In FMM, if particles
are located within each other’s near-field, their interaction is calculated directly. Otherwise,
a particle’s contribution to its nearest multipole is calculated and the interaction between
multipoles is then calculated by traversing up and down a tree.
A further example of an integral equation method that is of particular interest comes
from Alpert, Greengard and Hagstrom [5]. Here, the authors take the wave equation and
use a Laplace transform to create the equivalent PDE in the frequency domain. They then
take a centered finite difference of the temporal derivative followed by an inverse Laplace
transform to revert to the time domain. This gives that the finite difference approximation
of the temporal derivative is equivalent to the integral of the spatial derivative multiplied by
the appropriate Green’s function. This evolution scheme can be seen as an integral form of
the Lax-Wendroff scheme and the numerical properties of the scheme will depend upon the

6

technique used to compute the integral.
A less common approach is the Method of Lines Transpose (MOLT ), also known as the
transverse methods of lines or Rothe’s method or the horizontal line method, where the time
derivative is replaced by an algebraic approximation and the resulting system of Boundary
Value Problems (BVPs) must be solved. Here, each individual BVP represents the solution
at a given point in time. This approach has been relatively sparsely considered [8, 46, 56],
in part because the numerical solution to BVPs may be more challenging to find than for
the ODEs that result from the standard method of lines, especially in the case where the
BVP is stiff or contains boundary layers. While the MOLT may not be as widely used
as the MOL, it does not possess all of the same drawbacks and thus has been used in a
variety of applications such as solving the heat equation [21, 51], defect correction methods
[48], evolution equations [50], etc. With the continued development of boundary integral
methods [36, 55], the method of lines transpose approach may be very beneficial for certain
classes of problems. For example, Li et al. [55] constructed recurrence relations for efficiently
evaluating derivatives to the Green’s kernel for the non-oscillatory Helmholtz equation in
R3 .
In considering the vast array of previously developed temporal discretizations, two classes
of numerical schemes emerge: those which introduce numerical diffusion (i.e. dissipative
schemes), and those that do not (i.e. non-dissipative schemes). The non-dissipative schemes
will be useful for long time simulations that must maintain the strength of the wave amplitudes. The diffusive schemes are also of interest, for instance, in asymptotic preserving
applications in which the steady state approximation must be recovered.

7

1.2

Prerequisites

The following results that will be useful in showing the convergence of the MOLT scheme.

1.2.1

Eigenvalues of Perturbed Tri-diagonal Toeplitz Matrices

In section 3.1.2.1, we will use the inverse of a Toeplitz matrix to help determine stability
of the method. This inverse comes from Yueh and Cheng’s 2007 paper [70] on finding the
eigenvalues and inverses of a Tri-Diagonal Toeplitz matrix with perturbed corners. The
relevant matrix will be of the form



0
 b + γ a 0 ...


 a
b
a ...
0


 .
.
.. .. ..
.
An =  .
.
.
.
.
 .


 0
... a
b
a


0
... 0 a b + γ








.







Let λ be an eigenvalue of An with u its corresponding eigenvector. An u = λu can then
be written as

bu1 + au2 = λu1 − γu1
au1 + bu2 + au3 = λu2
.
.
.
aun−2 + bun−1 + aun = λun−1
aun−1 + bun = λun − γun

8

which is equivalent to

auk−1 + (b − λ)uk + auk+1 = 0,

k = 1, 2, . . . , n,

(1.3)

au0 = γu1 ,
aun+1 = γun .

Let ui , i ∈ {0, 1, . . . , n + 1} represent the first n + 1 terms of an infinite complex sequence
u = {ui }∞ . The equations in (1.3) then satisfy the recurrence relation
i=1

auk−1 + buk + auk+1 = λuk + fk ,

k∈N

under the condition

au0 = aun+1 = 0

with



−γu , k = 1


1




fk = −γun , k = n







0,

elsewise.
This relation is equivalent to

¯
¯ ¯
¯
ah2 + (b − λ)h + a u = (a¯1 + f )h
u

9

where

¯
h = {0, 1, 0, . . . },
¯
h2 = {0, 0, 1, 0, . . . },
a = {a, 0, 0, . . . },
¯
u1 = {u1 , 0, 0, . . . },
¯
f = {fk }∞ .
k=1

¯
¯ ¯
Note that h2 = h ∗ h, where ∗ is a convolution operator.
Since a = 0, this gives that

u=

¯
(a¯1 + f ) h
u
.
¯
¯ ¯
ah2 + (b − λ)h + a

(1.4)

¯
Consider the denominator as a quadratic in terms of h, whose roots are given by

−(b − λ) ±
δ± =
2a

√
ω

,

ω = (b − λ)2 − 4a2 .

Then δ+ δ− = 1 and thus,

δ± = cos(θ) ± i sin(θ)

for some θ in the strip {z ∈ C|0 ≤ Re(z) < π}.

10

(1.5)

Consider first the case where θ = 0. This will give that

cos(θ) =

λ−b
=⇒ λ = b + 2a cos(θ).
2a

Then, by partial fractions, equation (1.4) becomes

u=
=

1
ω

1
¯
δ− − h

−

1
¯
δ+ − h

¯
(a¯1 + f )h
u

2i
¯
{sin(jθ)} ∗ (a¯1 + f )h.
u
ω

Evaluating this convolution product will give

uj =

2i
n
{au1 sin(jθ) − δu1 sin((j − 1)θ) − Hj δun sin((j − n)θ)
ω

n
n
where Hj = 1 if j ≥ n and Hj = 0 if j < n.

This will then give

a sin(θ)un = au1 sin(nθ) − δu1 sin((n − 1)θ) and
0 = au1 sin((n + 1)θ) − δu1 sin(nθ) − δun sin(θ).

This can be rewritten as

Au1 + Bun = 0
Cu1 + Dun = 0.

11

(1.6)

Since u1 and un both can’t be zero, this implies that

A B
=0
C D
which gives the necessary condition

a2 sin((n + 1)θ) + δ 2 sin((n − 1)θ) − 2aδ sin(nθ) = 0.

For any θ that meets this qualification, the eigenvalues may then be found based on (1.6).
The case of θ = 0 will be similar, but now instead with equation (1.5) giving

1=

− (b − λ) ±

(b − λ)2 − 4a2
2a

=⇒ 2a = −(b − λ) ±

(b − λ)2 − 4a2

=⇒ (2a + (b − λ))2 = (b − λ)2 − 4a2
=⇒ 4a2 + 4a(b − λ) = −4a2
=⇒ b − λ = −2a =⇒ λ = b + 2a.

1.2.2

Schur Polynomial Test

The Schur polynomial test will be used in section 3.1.2.2 to guarantee that all the roots of
a polynomial are in the unit disc. A polynomial φd (z) =

d
i
i=0 ai z

is said to be a Schur

polynomial of exact degree d if φ(z) has exactly d roots, counting multiplicities, each of
which has magnitude less than one. The roots of this polynomial may be hard to find, but
the Schur polynomial test [67] may determine if all the roots are within the unit circle. To

12

check these conditions, define the corresponding polynomials
d

φ∗ (z)
d

ad−i z i and

=
i=0

φ∗ (0)φd (z) − φd (0)φ∗ (z)
d
.
φd−1 (z) = d
z

Theorem 1 (Schur Polynomial Test). φd (z) is a Schur polynomial of exact degree d if
φd−1 (z) is a Schur polynomial of exact degree d − 1 and |φd (0)| < |φ∗ (0)|.
d
Proof. Assume that φd−1 (z) is a Schur polynomial of exact degree d − 1 and that |φd (0)| <
|φ∗ (0)|.
d
First, note that on the unit circle itself,

|φ∗ (z)| = φd z −1
d

= |φd (z)|

since z −1 = z on the unit circle.
Define

ψ(z) =

zφd−1 (z)
.
φ∗ (0)
d

Then
φ∗ (0)φd (z) − φd (0)φ∗ (z)
d
|φd (z) − ψ(z)| = φd (z) − d
∗ (0)
φd
=

φd (0) ∗
∗
∗ (0) φd (z) < |φd (z)| = |φd (z)| on the unit circle.
φd

Thus, by Rouche’s theorm, φd (z) and ψ(z) must have the same number of zeros inside the
13

unit circle. Since φd−1 (z) is a Schur polynomial of exact degree d − 1, ψ(z) must have d
zeros, each of which lie within the unit circle, which implies that φd (z) also has d roots, each
of which lies within the unit circle.
Therefore, φd (z) is a Schur polynomial of exact degree d.

14

Chapter 2
MOLT Formulation and Integral Solution

2.1

Formulation

In the MOLT discretizations for the wave equation, the PDE is first discretized in time.
Both dissipative and non-dissipative classes of discretization schemes may be derived. The
following derivations are done in the context of one dimension but are easily extended to
higher dimensions.

2.1.1

Dissipative Schemes

To construct the dissipative schemes, first observe that the one-sided first and second order
n+1
finite difference approximations to the operator vtt = vtt (x, tn+1 ) are:

n+1
vtt =

v n+1 − 2v n + v n−1

+ O(∆t) and
∆t2
2v n+1 − 5v n + 4v n−1 − v n−2
n+1
vtt =
+ O ∆t2 .
2
∆t

15

(2.1)
(2.2)

Now let un = u(x, tn ) be the semi-discrete approximation to v(x, tn ), which satisfies the wave
equation (1.1). If the Laplacian term is evaluated implicitly, with either the substitution (2.1)
or (2.2), we obtain, respectively, the semi-discrete equations
un+1 − 2un + un−1
∆t2

= c2 un+1 and
xx

2un+1 − 5un + 4un−1 − un−2
∆t2

(2.3)

= c2 un+1 .
xx

(2.4)

Here and below, in accordance with context, we will reserve the variable v to mean the
continuous solution of the wave equation and u its semi-discrete approximation. It should
be emphasized that equations (2.3) and (2.4) produce numerical dissipation, which will be
shown in section 3.1. While not the focus of this work, this property is effective in recovering
magnetostatic limits and will be the subject of future work.

2.1.2

Purely Dispersive Schemes

A purely dispersive scheme may be constructed by carefully considering the choice of finite
n
difference approximations for both terms in vtt = c2 ∂xx v n , where the Laplacian term is

replaced with an appropriate time-averaged approximation. For example, since

∂xx v n = ∂xx

v n+1 + v n−1
2

+ O ∆t2 ,

a second order accurate semi-discrete approximation to equation (1.1) is
un+1 − 2un + un−1
(c∆t)2

= ∂xx

16

un+1 + un−1
2

.

(2.5)

As will be shown in section 3.1, this formulation removes numerical dissipation. This
purely dispersive scheme may be more favorable for long time simulations, which wish to
avoid the spurious effects of numerical diffusion.

2.2

Integral Solution

This work differs from other MOLT methods in that we consider integral solutions to the
boundary value problem, rather than utilizing finite difference or finite elements methods for
the spatial derivatives. Specifically, we solve the resulting BVP through a free space Green’s
function [46, 48, 56], which is then corrected with a boundary integral method.
Collecting un+1 in any of equations (2.3) - (2.5) results in

Lβj un+1 (x) = −

βj
(c∆t)2

fβj un , un−1 , un−2 , x ,

(2.6)

where the Helmholtz operator is

Lβj [u] :=

∂xx −

βj

u,

(c∆t)2

(2.7)

where βj is dependent upon the scheme chosen.
Equation (2.3) will correspond to choosing β1a = 1, for which the Helmholtz operator is

Lβ1a = ∂xx −

17

1
(c∆t)2

,

and

fβ1a = 2un (x) − un−1 (x).

Equation (2.4) will correspond to choosing β1b = 2, for which the Helmholtz operator is

Lβ

1b

= ∂xx −

2
(c∆t)2

and

fβ

1b

= 5un (x) − 4un−1 (x) + un−2 (x).

Similarly, equation (2.5) will also correspond to β2a = 2, so that the Helmholtz operator is
Lβ2a = Lβ

1b

and

fβ2a = 2un (x) −

(c∆t)2
L2 un−1 (x).
β2a

The subscripts will be dropped when discussing the general case.
Equation (2.6) can be solved formally by inverting the Helmholtz operator (2.7). The
resulting integral solution will consist of two components: a particular solution, up (x), which
solves the free space problem, and a homogeneous solution uh (x), which acts as a correction
to enforce the boundary conditions. This will produce a solution in the form

2
un+1 (x) = L−1 −αj fβj = un+1 (x) + un+1 (x),
p
β
h

18

(2.8)

where

αj =

βj
c∆t

,

and the respective particular and homogeneous solutions are given as

2
un+1 (x) = −αj
p

L
−L

fβj un (y), un−1 (y), un−2 (y) G(x|y) dy

(2.9)

and
un+1 (x) = un+1 (y)∂y G(x|y) − ∂y un+1 (y)G(x|y)
h

L
−L

.

(2.10)

The Green’s function for the Helmholtz operator (2.7) is




G(x|y) = − 1 e−α|x−y| ,
x ∈ R1


2α



1
2
G(x|y) = − 4 K0 (α|x − y|), x ∈ R





−α|x−y|


G(x|y) = e
,
x ∈ R3 ,
4π|x−y|

where K0 is a modified Bessel function.
Additionally, in the case where x ∈ R1 , the homogeneous solution (2.10) may be replaced
with the Ansatz
un+1 (x) = c1 e−α(L−x) + c2 e−α(L+x) ,
h
where c1 and c2 are determined by the imposed boundary conditions.

19

(2.11)

In the case of zero Dirichlet boundary conditions, this produces the system

un+1 (L) + c1 + c2 e−2αL = 0,
p
un+1 (−L) + c1 e−2αL + c2 = 0
p

which has a solution of
  
−1 

−2αL
n+1 (L)
e
c1   1
  −up

=
 
 

c2
e−2αL
1
−un+1 (−L)
p



−2αL
n+1 (L)
−e
1
 1
  −up

=


.
1 − e−4αL −e−2αL
n+1 (−L)
1
−up
Substituting in c1 and c2 in (2.11) then gives

un+1 (x) = −
h

sinh (α(L − x)) n+1
sinh (α(L + x)) n+1
up (−L) −
up (L).
sinh (2αL)
sinh (2αL)

In the case of zero Neumann boundary conditions, there is now the system

d n+1
u
(L) + αc1 − αc2 e−2αL = 0,
dx p
d n+1
up (−L) + αc1 e−2αL − αc2 = 0
dx

20

which has a solution of

  
−1 
d n+1
−αe−2αL   − up (L) 
c1   α

 =
  dx

 d
c2
αe−2αL
−α
− un+1 (−L)
dx p



d n+1
−2αL
−e
1
 1
  − dx up (L) 
.
=


 d

α(1 − e−4αL ) e−2αL
n+1 (−L)
−1
− up
dx
Note then that since u(x) is defined above in (2.9) as

un+1 (x) = −α2
p

L
−L

fβ un (y), un−1 (y), un−2 (y) G(x|y) dy.

This yields
L
d
d n+1
2
fβ un (y), un−1 (y), un−2 (y)
up = −α
G(x|y) dy
dx
dx
−L
L

= −α2
= −α2

−L
L
−L

fβ un (y), un−1 (y), un−2 (y) (−αsgn(x − y)G(x|y)) dy
−αsgn(x − y)fβ un (y), un−1 (y), un−2 (y) G(x|y) dy

and thus,
L
d n+1
up (L) = −α2
−αsgn(L − y)fβ un (y), un−1 (y), un−2 (y) G(L|y) dy
dx
−L

= −α3

L
−L

fβ un (y), un−1 (y), un−2 (y) G(L|y) dy

= −αun+1 (L)
p

21

while
L
d n+1
up (−L) = −α2
−αsgn(−L − y)fβ un (y), un−1 (y), un−2 (y) G(−L|y) dy
dx
−L
L

= α3

−L

fβ un (y), un−1 (y), un−2 (y) G(−L|y) dy

= αun+1 (−L).
p

This gives that
 
1
c1 
 =
1 − e−4αL
c2


 1

e−2αL

−e−2αL
−1



un+1 (L)
p





.

n+1 (−L)
−up

Substituting in for c1 and c2 in equation (2.11) gives this then as

un+1 (x) =
h

2.3

cosh (α(L − x)) n+1
cosh (α(L + x)) n+1
up (−L) +
up (L).
sinh (2αL)
sinh (2αL)

Numerical Approach

Consider first the discretization of the particular solution (2.9). Divide the domain into
non-overlapping regions Ωi , so that
N

un+1 (x) =
p
i=1

−α2

Ωi

fβ un (y), un−1 (y), un−2 (y) G(x|y) dy.

(2.12)

Generally, these regions may be of varying size, but in this dissertation, we will consider
a uniform partition for the composite midpoint rule, so that Ωi = xi−1/2 , xi+1/2 , i ∈
{1, 2, . . . N }, where xi = −L + (i − 1/2) ∆x and L = N ∆x/2. A midpoint quadrature

22

approximation to equation (2.12), evaluated at x = xj , gives for each j = 1, . . . N
N

un+1 (xj )
p

Aij fβ un (xi ), un−1 (xi ), un−2 (xi ) + EM (xj ),

=

(2.13)

i=1

where the remainder term satisfies
N

Bij ∂x fβ un (ξi ), un−1 (ξi ), un−2 (ξi )

EM (xj ) = ∆x
i=1

with ξi ∈ Ωi and

Aij = −2α2

xi+1/2
xi−1/2

xi+1/2

Bij = −2α2

G(xj |y) dy,

xi−1/2

(y − xi ) G(xj |y) dy.

Evaluating elements of the matrix A using the Green’s function gives


ν
 −|i−j|ν
e

sinh
, i=j
xi+1/2
−α|xj −y|
2
Aij = α
e
dy = 2

xi−1/2

1 − e−ν/2 ,

i=j
where
ν = α∆x =

β∆x
.
c∆t

Likewise, the elements of the matrix B satisfy

Bij = α

xi+1/2
xi−1/2

−α|xj −y|

(y − xi ) e

dy =





c(ν)e−|i−j|ν , i = j


0,


23

i=j

(2.14)

where c(ν) = cosh

ν
2
ν
− sinh
. In practice, we are developing implicit methods, and
2
ν
2

as such, the case where ∆x < c∆t, corresponding to ν < 1, is of particular interest. Taking
a Taylor expansion for small ν shows

2
ν
ν
− sinh
c(ν) = cosh
2
ν
2

∞

=
n=0

=1+
=

(ν/2)2n 2
−
(2n)!
ν

∞
n=0
4

(ν/2)2n+1
(2n + 1)!

(ν/2)2 (ν/2)
2
+
−
2
24
ν

ν (ν/2)3 (ν/2)5
+
+
2
6
120

+ O ν6

ν2
ν4
+
+ O ν6 .
12 480

Thus, the error is bounded by
∞

EM (xj ) ∞ ≤ 2∆xc(ν) ∂x fβ ∞

e−iν

(2.15)

i=1
2e−ν

= ∆xc(ν) ∂x fβ ∞
1 − e−ν
ν∆x
∂x fβ ∞ .
≤
6

Here, we must note that each term in fβ is a numerical solution to the wave equation at a
previous time step. Thus, since

∂
1∂
f =±
f and we have already discretized the temporal
∂x
c ∂t

derivative, we have that ∂x fβ behaves like fβ ∆t−1 , giving that the error from equation (2.15)
is actually O ν 2 .
For a trapezoidal approximation, consider the points xi = −L+i∆x with Ωi = [xi−1 , xi ].
This approximation then gives for each j = 0, 1, . . . N
N

un+1 (xj ) =
p

Cij fβ un (xi ), un−1 (xi ), un−2 (xi ) + ET (xj ),
i=0

24

(2.16)

where the error term is
∆x2
ET (xj ) =
2

N

Dij ∂xx fβ un (ξi ), un−1 (ξi ), un−2 (ξi ) .
i=0

Analogous evaluation of the matrix terms gives

Cij = α

=

xi
xi−1

y − xi−1
∆x

−α|xj −y|
dy

e

xi+1

+α

xi+1 − y
∆x

xi

e

−α|xj −y|





2
−|i−j|ν , i = j
2 2 sinh (ν/2)e
ν
 −ν
e − 1 + ν,


dy

(2.17)

i=j

and

Dij = α

xi
xi−1

where again c(ν) = cosh

(y − xi−1 ) (y − xi ) e

−α|xj −y|

2
ν

dy =

c(ν)e−|i−j+1/2|ν

2
ν
ν
− sinh
. The corresponding error bound is
2
ν
2

ET (xj ) ≤
=
≤
≤

∆x2
2
∆x2
ν
∆x2
ν

∂xx fβ ∞

∞

2
ν

e−|i−j+1/2|ν

c(ν)

(2.18)

i=−∞
∞

c(ν) ∂xx fβ ∞ 2

e−|i−j|ν e−ν/2

i=0
2e−ν/2

c(ν) ∂xx fβ ∞
1 − e−ν

∆x2
∂xx fβ ∞ .
3

Similar to the proof of the error estimate of the midpoint quadrature, ∂xx fβ behaves like
fβ ∆t−2 , showing that the error in equation (2.18) is truly O ν 2 .
25

Since both of the error estimates for the midpoint and trapezoidal quadrature rule are
O ν 2 , this will demand that we must satisfy that our spatial step size goes to zero at an
exponentially faster rate than the temporal step size. As we are creating an implicit scheme,
this is much less of a restriction than for an explicit scheme. However, in section 3.1.1.3,
we will show that the method of undetermined coefficients may be used to find weights that
will change the error from O ν 2 to O ∆x2 , as is desired.

26

Chapter 3
Analysis of the Free Space Solution

3.1

Convergence

The main result of this chapter is the following convergence theorem.
Theorem 2. Let un (x) be the semi-discrete MOLT solution satisfying equations (2.8), (2.9),
(2.10), with initial and boundary conditions specified by equation (1.1). Then the semidiscrete solutions satisfying equations (2.3), (2.4) and (2.5) will converge to v(x, tn ) as
O ∆tβ , with β = 1, 2, 2 respectively. Further, suppose the fully discrete solution un is
j
computed using either the midpoint (2.13) or trapezoidal (2.16) quadrature rules. Then, un
j
will converge to v(xj , tn ) as O ∆tβ + ν 2 and will be stable for 0 ≤ t ≤ T .
The proof of this theorem follows from the Lax–Richtmyer equivalence theorem upon
establishing consistency and stability. Consistency is established in section 3.1.1 and stability
in 3.1.2. In section 3.1.1.3, we will determine quadrature weights that will create a fully
discrete scheme that is O ∆tβ + ∆x2

In section 3.2, the dispersion of these methods is

discussed.
27

3.1.1

Consistency of the MOLT method

3.1.1.1

Non-dissipative Scheme

Let v be the the continuous solution satisfying the wave equation (1.1) with L = ∞. Now,
the integral solution for the semi-discrete solution for (2.5) in free space is given by

un+1 (x) = I [un ] (x) − un−1 (x),

where I [v n ] is defined by

I [v n ]

∞

=

2α2 G(x|y)v n (y)dy

∞

=α

e−α|z| v n (x + z) dz

−∞

−∞

after applying the change of variables y − x = z.
Inserting v into this equation, we can define the global truncation error as

τn =

1
(c∆t)2

v n+1 + v n−1 − I [v n ] ,

The integral can be split at z = 0, giving
∞

I = (I+ ) + (I− ) = α

0

e−αz v n (x + z)dz

28

∞

+α

0

e−αz v n (x − z)dz.

(3.1)

For each integral, repeatedly apply integration by parts to get
∞

I± = α

0

e−αz v n (x ± z) dz
∞

∞

= −(±1)v n (x ± z)e−αz 0 +

0

±∂x v n (x ± z)e−αz dz

.
.
.


k

= −α
j=0
∞

+α

0

∞

±1
∂x
α

j

±1
∂x
α

k+1

v n (x ± z)e−αz 
0

v n (x ± z)e−αz dz.

Notice that when evaluating at z = 0, the alternating signs cause the odd derivatives to
cancel and the even derivatives to combine for I = I+ + I− . Furthermore, the upper limit
makes no contribution due to the dominance of the negative exponential.
Truncating the integration by parts at the leading order for the error, corresponding here
to k = 3, gives

I [v n ]

∞
1
1
1
n (x + z)e−αz
∂xx +
∂xxx v
=
1 + ∂x +
α
α2
α3
0
∞
1
1
1
∂xx −
∂xxx v n (x − z)e−αz
+ 1 − ∂x +
α
α2
α3
0
∞

+α

0

= 2 vn +

1
α4

n
n
(vxxxx (x + z) + vxxxx (x − z)) e−αz dz

1 n
vxx
α2

∞

+α

n
= 2v n + (c∆t)2 vxx +

0

1
α4

n
n
(vxxxx (x + z) + vxxxx (x − z)) e−αz dz

∞
(c∆t)4
n
(2vxxxx (ξ)) α
e−αz dz.
4
0

(3.2)

Note that we use the mean value theorem for integrals in the last line. The truncation

29

error (3.1) then becomes

τn =
=
=

1
(c∆t)2
1
c2
1

n
v n+1 + v n−1 − 2v n − (c∆t)2 vxx −

n
vtt +

∆t2 n
v + O ∆t4
12 tttt

∆t2 n
vtttt + O ∆t4
12

−

c2
−5(c∆t)2 n
=
vxxxx (ζ) + O ∆t4
12

n
− vxx −

(c∆t)4 n
vxxxx (ξ)
2

(c∆t)2 n
vxxxx (ξ)
2

(c∆t)2 n
vxxxx (ξ)
2

where the terms v n±1 have been expanded about tn and differentiated in the context of the
wave equation to determine that

vtttt = c2 vttxx = c4 vxxxx .

Hence, the semi-discrete error is second order in time.

3.1.1.2

Dissipative Schemes

The same ideas may be applied to the dissipative free space solver (2.6). However, the
expansion is now about tn+1 , where now the first order semi-discrete solution is given by

1
un+1 (x) = I [un ] (x) + I un−1 (x)
2

30

where I [u] is defined, as before, by
∞

I [u] =

2α2 G(x|y)u(y)dy = α

−∞

∞

e−α|z| u(x + z) dz.

−∞

Substituting in the continuous solution v and once again changing variables to z = y − x,
the global truncation error is defined as

τ n+1 =

1
(c∆t)2

1
v n+1 − I [v n ] − I v n−1
2

.

Repeating the integral splitting and integration by parts done just above, it becomes easy
to see that truncating after three steps gives

1
I [v n ] − I v n−1
2
∞
1
1
1
1
=
1 + ∂x +
∂xx +
∂xxx
v n − v n−1 (x + z)e−αz
α
2
α2
α3
0
∞
1
1
1
1
∂xx −
∂xxx
+ 1 − ∂x +
v n − v n−1 (x − z)e−αz
α
2
α2
α3
0
∞

=2

1

∞

+α

1

1
1
(x + z) + v n − v n−1
(x − z) e−αz dz
v n − v n−1
4
2
2
0 α
xxxx
xxxx
1
1
1
v n − v n−1 +
v n − v n−1
2
2
2
α
xx

1
1
v n − v n−1
(x + z) + v n − v n−1
(x − z) e−αz dz
4
2
2
0 α
xxxx
xxxx
1
1
= 2 v n − v n−1 + (c∆t)2 v n − v n−1
2
2
xx
+α

+

(c∆t)4
4

∞
1
2 v n − v n−1
(ξ) α
e−αz dz.
2
0
xxxx

31

Thus, the global truncation error is

τ n+1 =
=

1
(c∆t)2
1

xx

+ O ∆t4

n+1
n+1 − v n+1 + O ∆t2
vtt − ∆tvttt
xx

c2

=−

v n+1 − 2v n + v n−1 − (c∆t)2 2v n − v n−1

1
c2

n+1
∆tvttt (ζ) + O ∆t2

n+1
= −∆tvtxx (ζ) + O ∆t2 ,

which shows that this method is first order in time.
The second order semi-discrete solution is given by

1
5
un+1 (x) = I [un ](x) + I un−1 (x) − I un−2 (x),
4
4

and thus

τ n+1 =

1
(c∆t)2

5
1
v n+1 − I [v n ] + I v n−1 − I v n−2
4
4

.

Using the same calculations as above,

I [v n ] = 2 v n +
=

2v n

1 n
vxx
α2

n
+ (c∆t)2 vxx

∞

+α

0

1
α4

n
n
(vxxxx (x + z) + vxxxx (x − z))e−αz dz

∞
(c∆t)4
n
+
(2vxxxx (ξ)) α
e−αz dz.
4
0

32

and therefore, after utilizing the linearity of I, this shows that

τ n+1 =

1

5
5
n
n−1
v n+1 − v n − (c∆t)2 vxx + 2v n−1 + (c∆t)2 vxx
2
2
4
(c∆t)

1 n−2
1
− v n−2 − (c∆t)2 vxx + O ∆t4
2
4
=

1

1 n+1
1 n+1 3∆t2 n+1
vtt
− vxx −
vxxtt + O ∆t4
2
2
4

c2
3(c∆t)2 n+1
=−
vxxxx + O ∆t4 ,
4

which gives that the method is indeed second order.

3.1.1.3

Fully Discrete Schemes

In the case of the fully discrete implementations, the results from sections 3.1.1.1 and 3.1.1.2
will still follow in the exact same manner, with the addition that I [u] now contributes an
O ν 2 error caused by the quadrature estimate.
A quick observation from equations (2.14) and (2.17) will show that the terms that come
from integrating equation (2.12) using quadrature will be of the form
∞

e−|k−j|ν u(xj−k ) + u(xj−k ) .

I [u](xj ) = w0 u(xj ) + w1
k=1

Consider now the truncation error given by the non-dissipative scheme, shown in (3.1),
adapted for the fully discrete scheme as

τ n (xj ) =

1
(c∆t)2

un+1 (xj ) + un−1 (xj ) − I [un ](xj ) .

33

(3.3)

Let

2w
w0 = 2 − ν 1 = 2 − w1
e −1

∞

2e−kν
k=1

(3.3) then becomes

τ n (x

j)

=

un+1 xj − 2un xj + un−1 xj
∆t2

= un (xj ) +
tt
= un (xj ) +
tt

∆t2
12

utttt (ξ1 ) −

w1
∆t2

∆t2
utttt (ξ1 ) − w1
12

w
− 12
∆t

∞

k=1

e−kν un (xj+k ) − 2un (xj ) + un (xj−k )
k=1

(k∆x)2 un (xj ) +
xx

e−kν

k=1
∞
e−kν

∞

(k∆x)4 n
uxxxx (ξ2 )
12

(kν)4 ∆t2 n
(kν)2 n
utt (xj ) +
utttt (ξ2 )
2
48

= η1 un (xj ) + ∆t2 η2 un (ξ)
tt
tttt

(3.4)

where
w ν2
η1 = 1 − 1
2
η2 =

∞

e−kν k 2 ,

k=1
4 ∞
ν

w
1
− 1
12
48

e−kν k 4 .

k=1

Note that
∞

e−kν k j = (−1)j
k=1

∞
d j
e−kν = (−1)j
dν
k=1

34

d j 1
,
dν
eν − 1

most notably with
∞
k=1
∞
k=1

1
ν
ν
e−kν k 2 = coth
csch2
4
2
2

and

1
ν
ν
e−kν k 4 = coth
csch2
4
2
2

3csch2

ν
+1 .
2

We must choose w1 such that η1 = 0 if we wish to have the order be O ∆x2 . Thus,
ν
ν
w1 ν 2 1
coth
csch2
=0
2 4
2
2
8
ν
ν
=⇒ w1 = 2 tanh
sinh2
.
2
2
ν

η1 = 0 =⇒ 1 −

Then,
8ν 4
ν
ν 1
ν
ν
1
−
tanh
sinh2
coth
csch2
12
48
2
2 4
2
2
2
1
ν
ν
=
−
3csch2
+1 .
12 24
2

η2 =

3csch2

Using the expansion

csch2 (x)

1
1
= 2− +
3
x

1
+
36

7 2
x2 + O x4
360

gives that

η2 =

−5 3ν 2
−
+ O ν4
12
24

35

ν
+1
2

and thus, substituting in η2 into (3.4) gives that the truncation error is now

τ n = ∆t2

−5 3ν 2
−
+ O ν4
12
24

= O ∆t2 + ∆x2 .

Similar results will hold for the same w0 and w1 for the dissipative schemes.

3.1.2

Stability of the MOLT method

To establish the stability of the MOLT method, we treat the fully discretized system as a
one step method. To that extent, define

un = [un (x1 ), . . . , un (xN )]T ,
¯

n ≥ 1,

where xj = −L + (j − 0.5) ∆x in the case of the midpoint rule and xj = −L + j∆x when
considering the trapezoidal rule. Then, for the dissipative scheme corresponding to equation
(2.3),

1
¯
un+1 = At un − un−1 .
¯
¯
2

(3.5)

For the second order dissipative scheme in accordance with equation (2.4), we have

1
u
u
¯
un+1 = At 5¯n − 4¯n−1 + un−2
¯
4

(3.6)

and for the non-dissipative scheme, corresponding to equation (2.5),

un+1 = At un − un−1 .
¯
¯
¯
36

(3.7)

In each of the cases, the matrix At will possess the same structure, although the individual
matrix elements, which are based upon ν, will differ. At will consist of a particular and
homogeneous contribution, in analog to the free space and boundary terms in the integral
solution of the continuous problem. It also contains the information of the discretization
scheme in use; in this case, the midpoint or trapezoidal quadrature rules. If we write At =
Ap + Ah , then for the midpoint rule, Ap is given by the matrix A in equation (2.14), and
the matrix C of equation (2.17) in the case of the trapezoidal rule. The boundary matrix
entries are

(Ah )ij = −2a0 sinh
+

ν
2

sinh (α(L + xi )) −α L−xj
e
sinh (2αL)

sinh (α(L − xi )) −α L+xj
e
sinh (2αL)

(3.8)

,

where the coefficient a0 distinguishes between the midpoint and trapezoidal rules with

a0 =




1,


Midpoint Rule
(3.9)

2

 sinh ν , Trapezoidal Rule.

ν
2

We will establish stability for these schemes by studying the eigenvalues of Ap and Ah .
Since Ah will be exponentially small away from the boundaries, the main part of this work
will focus on analyzing the eigenvalues of the free-space problem, which would neglect Ah
altogether. Once the free space problem is understood, we will proceed with the boundary
correction matrix.

37

3.1.2.1

The Eigenvalues of Ap

In order to simplify the proceeding analysis, which is for two different methods of discretization, write the matrices Ap = A of equation (2.14) and Ap = C of equation (2.17) in the
form
ν
T,
2

Ap = 2˜0 IN + 2˜1 sinh
a
a

(3.10)

where IN is the N × N identity matrix, and the constants a0 and a1 distinguish the quadra˜
˜
ture rule. They are given by



1 − e−ν/2 − sinh ν ,

2
a0 , a1 =
˜ ˜
1

ν


e−ν − 1 + ν − 2 sinh2
ν
2

1,

Midpoint Rule
(3.11)

,

2
ν
sinh
, Trapezoidal Rule.
ν
2

Note that these expressions can be simplified to 1−cosh

ν
1
and 1− sinh (ν), respectively.
2
ν

Finally, the matrix T is a Toeplitz matrix and is given by

Tij = x|i−j| ,

x = e−ν .

Notice that the eigenvalues of Ap can be found from the eigenvalues tk of T , since Ap differs
from T by a scaling factor and a multiple of the identity matrix. Therefore, the eigenvalues
λk of Ap are given by

λk = 2˜0 + 2˜1 sinh
a
a

38

ν
t .
2 k

(3.12)

Now, due to its structural symmetry, the inverse matrix T −1 has a special form, and the
eigenvalues may be found analytically. In fact, T −1 is a nearly-Toeplitz, tri-diagonal matrix,






1 

−1 =
T

2
1−x 






1

−x

0

...

0 


2 −x
−x 1 + x
...
0 


.
. ,
...
...
..
.
. 
.
.
. 

2 −x 

0
...
−x 1 + x


0
...
0
−x
1

x = e−ν .

Since this is a tri-diagonal Toeplitz matrix with permuted corners, recall that based on the
results in section 1.2.1, the eigenvalues of this matrix are given by

t−1 =
k

x2 − 2x cos (θk ) + 1
1 − x2

,

where θk is a root of

sin((N + 1)θ) + x2 sin((N − 1)θ) − xsin(N θ) = 0,

θ ∈ (0, π).

This gives that the eigenvalues of our Toeplitz matrix are given by

tk =

1 − x2
x2 − 2x cos (θk ) + 1

39

,

Making use of x = e−ν , we can find an explicit form for the eigenvalues by solving

tk =
=
=

1 − x2
x2 − 2x cos (θk ) + 1
1 − e−2ν
e−2ν − 2e−ν cos(θk ) + 1
eν − e−ν
e−ν − 2cos(θk ) + eν

=

2sinh(ν)
2cosh(ν) − 2cos(θk )

=

sinh(ν)
cosh(ν) − cos(θk )

=
=

2sinh (ν/2) cosh (ν/2)
1 + 2sinh2 (ν/2) − 1 + 2sin2 (θk /2)
sinh (ν/2) cosh (ν/2)
sinh2 (ν/2) + sin2 (θk /2)

.

Thus, the eigenvalues are given by

tk = cosh (ν/2)

sinh (ν/2)
sinh2 (ν/2) + sin2 (θk /2)

,

0 < θk < π.

This leads to the following lemma.
Lemma 3. The eigenvalues of the matrix given by equations (2.14) and (2.17) are contained
in the interval (0, 2).
Proof. This is a consequence of equation (3.12), using the coefficients from (3.11). In the

40

former case, the eigenvalues for k = 1, 2, . . . N , are
sinh2 (ν/2)

λk = 2 1 − cosh (ν/2) + cosh (ν/2)
= 2 1 − cosh (ν/2)

sinh2 (ν/2) + sin2 (θk /2)

sin2 (θk /2)

,

sinh2 (ν/2) + sin2 (θk /2)

0 < θk < π.

(3.13)

We obtain the bounds by letting θk = 0, π, and therefore

2 − 2sech

ν
2

< λk < 2,

(3.14)

where θk = 0 corresponds to 2 and θk = π corresponds with 2 − 2sech

ν
. Similarly, for
2

the trapezoidal rule, the eigenvalues are for k = 1, 2, . . . N ,

λk = 2 1 −
=2 1−

1
2
sinh (ν) + sinh (ν/2) cosh (ν/2)
ν
ν
1
sinh (ν)
ν

sinh2 (ν/2)
sinh2 (ν/2) + sin2 (θk /2)

sin2 (θk /2)

sinh2 (ν/2) + sin2 (θk /2)

,

0 < θk < π.

(3.15)

Now, the bounds obtained by letting θk = 0, π, are

2−

4
ν
tanh
ν
2

< λk < 2,

where θk = 0 corresponds to 2 and θk = π corresponds with 2 −

(3.16)

4
ν
tanh
.
ν
2

The lower bound in (3.14) and (3.16) will be identically zero only for ν = 0 and positive
otherwise. Therefore, the lemma is proven.
Lemma 4. The matrices A and C given by equation (2.14) and (2.17) respectively, are
41

symmetric positive definite.
Proof. Note first that Ap is symmetric since T is symmetric as Tij = x|i−j| = x|j−i| = Tji .
Also,

sinh(2αL)
(A ) =
−2a0 sinh (ν/2) h ij
−α(L−xj )

sinh(α(L + xi ))e

+ sinh(α(L − xi ))e

= eα(L+xi ) − e−α(L+xi ) e
α(xi +xj )

=e

= e

− e−2αL e

α(L+xj )

−e

−α(L−xj )

α(−xi +xj )

−α(L+xj )

−α(L+xj )

+ eα(L−xi ) − e−α(L−xi ) e

−α(xi +xj )

+e

− e−2αL e

α(L−xj )

e−α(L−xi ) + e

−α(L+xj )

α(xi −xj )

−α(L−xj )

−e

e−α(L+xi )

= sinh(α(L + xj ))e−α(L−xi ) + sinh(α(L − xj ))e−α(L+xi )
=

sinh(2αL)
(A ) .
−2a0 sinh (ν/2) h ji

Thus, the matrices A and C are symmetric. They are positive definite because of the previous
lemma which proved that all of the eigenvalues were positive.

3.1.2.2

Dissipative Schemes

We are now in a position to show stability for the free-space problem using Von Neumann
analysis. Let the numerical solution be given by un = µn u0 for n ≥ 0. Inserting this
¯
¯
expression into equations (3.5) and (3.6) with At = Ap , and canceling the common factor
produces the stability equations

µ2 u0 =
¯
µ3 u0 =
¯

µ−

1
2

Ap u0 and
¯

1
5µ2 − 4µ + 1 Ap u0 .
¯
4
42

(3.17)

Now, since Ap is a symmetric positive definite matrix, it is diagonalizable as Ap = QΛQ−1 ,
where Q is orthonormal and Λ is the diagonal matrix with entries given by (3.13) or (3.15).
Thus, solvability of the expression (3.17) leads to the degree 2N Von-Neumann polynomial

p(µ) = µ2 IN − µ −

1
2

Ap

= µ2 QIN Q−1 − µ −
1
= µ2 IN − µ −
2

1
2

QΛQ−1

Λ = 0.

Since this latter form is expressed as the determinant of a diagonal matrix, the Von-Neumann
polynomial can be written in factored form, and the amplification factors µ± will be deterk
mined in pairs by the expression

λ
µ2 − λk µk + k = 0,
k
2

k = 1, 2, . . . N.

Solving these resulting quadratic equations gives

µ±
k

=

λk ±

λ2 − 2λk
k
2

.

Since λk ∈ (0, 2), it is immediately apparent that the amplification factors form complex
conjugate pairs. Further, since each θk will be distinct, this means that for each k, the roots
of the Von-Neumann polynomial will also be distinct. If it is further shown that |µk | ≤ 1,
then stability of the numerical scheme will be established. Taking the modulus of the roots,

43

we have

µ± =
k

λk 2 λk
−
+
2
2

λk 2
=
2

λk
< 1.
2

Thus, the first order system will be stable, and, since |µ| < 1, dissipative.
Repeating the process for the second order formulation gives

4µ3 − 5λk µ2 + 4λk µk − λk = 0.
k
k

Lemma 5. The magnitude of the roots of 4µ3 − 5λk µ2 + 4λk µk − λk are all less than one.
k
k
Proof. This will be shown using the Schur polynomial test given in section 1.2.2.
Define

φ0 (z) = 4z 3 − 5λk z 2 + 4λk z − λk and
φ∗ (z) = −λk z 3 + 4λk z 2 − λk z + 4.
0

Note that |φ0 (0)| = λk < 4 = |φ∗ (0)|.
0
Then,
4 4z 3 − 5λk z 2 + 4λk z − λk + λk −λk z 3 + 4λk z 2 − λk z + 4
φ1 (z) =
z
= (16 − λk ) z 2 + 4λ2 − 20λk z + 16λk − 5λ2 and thus
k
k
φ∗ (z) = 16λk − 5λ2 z 2 + 4λ2 − 20λk z + 16 − λk .
1
k
k

This gives that φ1 (0) = 16λk − 5λ2 > 0 and φ∗ (0) = 16 − λk > 0. Thus, we can check the
1
k

44

condition of |φ1 (0) ≤ |φ∗ (0)| by finding the roots of φm (0) − φ∗ (0) = −5λ2 + 17λk − 16.
m
1
k
Solving this quadratic shows that the roots are imaginary and thus the equation is negative
everywhere since using λk = 0 gives −16.
Taking this further by one step yields

φ2 (z) =

1
z

(16 − λk ) (16 − λk ) z 2 + 4λ2 − 20λk z + 16λk − 5λ2
k
k

− 16λk − 5λ2
k
=

16λk − 5λ2 z 2 + 4λ2 − 20λk z + (16 − λk )
k
k

(16 − λk )2 − 16λk − 5λ2
k

2

z + 16 − 17λk + 5λ2
k

4λ2 − 20λk ,
k

which is obviously a polynomial of degree 1.
Furthermore, this yields that

˜
φ2 (z) =

(16 − λk )2 − 16λk − 5λ2
k

16 − 17λk + 5λ2
k

2

4λ2 − 20λk
k

=

4λk 5λ3 − 42λ2 + 85λk − 80
k
k
16 − 17λk + 5λ2
k

4λ2 − 20λk
k

.

˜
Mathematica can then be used to show that φ2 (z) has a critical point at 0 and no other
˜
critical points in [0, 2]. Thus, since |φ(1)| ≈ 0.4848, I have that φ2 is a Schur polynomial
and thus subsequently, so are φ1 and φ0 .

3.1.2.3

Non-dissipative Scheme

The process of the previous section is now repeated on (3.7). This gives rise to a polynomial
of the form

µ2 − λk µk + 1 = 0,
k
45

k = 1, 2, . . . N.

Thus, the amplification factors are now given as

µ±
k

=

λ2 − 4
k

λk ±
2

.

This time, the modulus of each complex conjugate pair is

|µ± |
k

=

λk 2
+1−
2

λk 2
= 1.
2

Thus, since the eigenvalues are all distinct, the method is stable, and, since |µ| = 1, the
scheme is non-dissipative.

3.2

Dispersion

We now analyze the phase error of the free-space approximation. The continuous dispersion
relation results from looking at sinusoidal solutions of the wave equation

v(x, t) = ei(kx−ωt) with ω 2 = c2 k 2 .

We now analyze the semi-discrete dispersion relations and define the phase error for the free
space dissipative and non-dissipative schemes.

46

3.2.1

Semi-discrete Non-dissipative Scheme

ω
For the semi-discrete equation, define un = ueikx−i˜ (n∆t) , where ω denotes the discrete
ˆ
˜

temporal frequency. Substituting this ansatz into the non-dissipative scheme yields
∞

ω
ω
ω
ˆ
e−i˜ ∆t + ei˜ ∆t ueikx−i˜ (n∆t) = α

ω
e−α|x−y| ueiky−i˜ (n∆t) dy.
ˆ

−∞

Cancel the common term ue−iωn∆t to get
ˆ
α ∞ −α|x−y|−ik(x−y)
e
dy
2 −∞
α ∞ −αz ikz
e
e + e−ikz dz
=
2 0

cos (˜ ∆t) =
ω

(3.18)

after the change of variables z = y − x. Evaluation of the integral gives

cos (˜ ∆t) =
ω

α2
α2 + k 2

=

1
(kc∆t)2
1+
2

.

(3.19)

This is the semi-discrete analog to the continuous dispersion relation. To avoid aliasing,
consider only the wave numbers in the region 0 ≤ ω ∆t ≤ π. The phase error can then be
˜
defined as

Φ(∆t) =

ω
˜
−1 .
kc

47

For a fixed frequency k, we can analyze the convergence of this approximation using a Taylor
series approximation for small kc∆t. First, make use of the trigonometric identity

cos (2θ) = 1 − 2 sin2 (θ)

in (3.19) to obtain

1 − 2 sin2

=⇒ 2 sin2

=⇒ 2 sin2

ω∆t
2

=

ω∆t
2

=1−

1
(kc∆t)2
2
1

1+

(kc∆t)2
2
2
(kc∆t)
=
,
2 + (kc∆t)2
1+

ω∆t
2

which gives

sin2

where z =

ω ∆t
˜
2

=

(kc∆t)2 


4


2

= z
,
(kc∆t)2  1 + 2z 2
1+
2

1

kc∆t
. Taking the square root and arcsin of both sides of the equation now yields
2

ω=
˜

2
arcsin
∆t

48

z
1 + 2z 2

.

(3.20)

The first terms of the binomial expansion for the denominator are given by

1 + 2z 2

−1/2

2
1
3
= 1 − 2z 2 +
2z 2
2
8
2
(kc∆t)
3(kc∆t)4
=1−
+
4
32

which, in turn, gives that

z
1 + 2z 2

=

kc∆t
2

1−

(kc∆t)2
4

+ O ∆t5 .

Substituting this into the arcsin expansion

arcsin(z) = z +

z 3 3z 5
+
+ O z7
6
40

now turns (3.20) into

ω=
˜

2  kc∆t


∆t  2

= kc 1 −

1−

(kc∆t)2
4

+

5
(kc∆t)2 + O ∆t4
24

kc∆t
2

(kc∆t)2
1−
4
6

,

from which it follows that the phase error is second order.

49

3




+ O ∆t5 


3.2.2

Semi-discrete Dissipative Schemes

ω
For the dissipative scheme, similarly define un = ueikx−i˜ (n∆t) . Substituting this ansatz
ˆ

into the first order dissipative scheme yields

ω
e−i˜ (n+1)∆t ueikx =
ˆ

α ∞ −i˜ n∆t iky
ω
e ω
ue − e−i˜ (n−1)∆t ueiky e−α|x−y| dy.
ˆ
ˆ
2 −∞

Once again, cancel the common term ue−i˜ n∆t to get
ˆ ω
α ∞
ω ˆ
ueik(y−x) − ei˜ ∆t ueik(y−x) e−α|x−y| dy
ˆ
2 −∞
α ∞ ik(y−x)−α|x−y|
ω
e
2 − ei˜ ∆t dy
=
2 −∞

ω
e−i˜ ∆t =

=⇒

ω
e−i˜ ∆t
ω
2 − ei˜ ∆t

=

α ∞ ik(y−x)−α|x−y|
e
dy.
2 −∞

Use the change of variables z = y − x to get
ω
e−i˜ ∆t
ω
2 − ei˜ ∆t

=

α ∞ −αz ikz
e
e + e−ikz
2 0

dz.

Performing the integration gives
ω
e−i˜ ∆t
ω
2 − ei˜ ∆t

=

ω
=⇒ e−i˜ ∆t =

α2
α2 + k 2
ω
2 − ei˜ ∆t
1 + (kc∆t)2

.

From this expression we would like to isolate ω , and so we solve the corresponding
˜

50

ω
quadratic polynomial for λ = ei˜ ∆t ,

λ2 − 2λ + 1 + (kc∆t)2 = 0,

which admits complex solutions λ = 1 + ±ikc∆t. Thus, the temporal frequencies will be
complex (dissipative), and are of the form

1
log (1 + ikc∆t)
i∆t
ikc∆t
= kc 1 −
+ O ∆t2
2

ω=
˜

,

which shows that the phase error is first order.

ω
We then repeat this for the final method by putting the ansatz un = ueikx−i˜ (n∆t) into
ˆ

the second order dissipative scheme. This gives

ω
e−i˜ (n+1)∆t ueikx =
ˆ

α ∞
2 −∞

1 ω
5 −i˜ n∆t iky
ω
e ω
ue − 2e−i˜ (n−1)∆t ueiky + e−i˜ (n−2)∆t ueiky e−α|x−y| dy.
ˆ
ˆ
ˆ
2
2

Cancel the common term ue−i˜ n∆t to get
ˆ ω
α ∞
2 −∞
α ∞
=
2 −∞

ω
e−i˜ ∆t =

=⇒

5 ik(y−x)
1 ω
ω
e
− 2ei˜ ∆t eik(y−x) + e2i˜ ∆t eik(y−x) e−α|x−y| dy
2
2
5
1 ω
ω
− 2ei˜ ∆t + e2i˜ ∆t eik(y−x)−α|x−y| dy
2
2

ω
2e−i˜ ∆t
ω
5 − 4ei˜ ∆t

ω
+ e2i˜ ∆t

∞

=

eik(y−x)−α|x−y| dy.

−∞

51

Again, use the change of variables z = y − x to get
ω
2e−i˜ ∆t
ω
ω
5 − 4ei˜ ∆t + e2i˜ ∆t

=

α ∞ −αz ikz
e
e + e−ikz
2 0

dz.

Performing the integration gives
ω
2e−i˜ ∆t
ω
ω
5 − 4ei˜ ∆t + e2i˜ ∆t

=

ω
=⇒ 2e−i˜ ∆t =

α2
α2 + k 2
ω
ω
5 − 4ei˜ ∆t + e2i˜ ∆t
(kc∆t)2
1+
2

.

ω
Setting λ = ei˜ ∆t and setting this as a polynomial gives

λ3 − 4λ2 + 5λ − 2 − (kc∆t)2 = 0.

Solving this in Mathematica gives the solutions as

λ=

√
√
√
√
√
3
2
4
R 4 1 − i 3 (1 + i 3)R 4 1 + i 3 (1 − i 3)R
+
+ √ , − √
−
, − √
−
√
√
3
3R
3
332 3
3 3 4R
632
3 3 4R
632

where R =

3

, (3.21)

√
2 + 27(kc∆t)2 + 3 3 4(kc∆t)2 + 27(kc∆t)4 .

From here, consider the asymptotic limit as ∆t → 0, which gives that these roots approach
√
3
4
2
R
respectively {2, 1, 1}. Note that the first root, +
+ √ is real and thus corresponds
3
3R
332
to a parasitic mode instead of a propagating mode and thus it will not be considered any

52

further. For the two complex roots, consider perturbative roots of the form

λ = 1 + δ,

δ = c1 + c2 2 + c3 3 + . . . ,

= kc∆t.

(3.22)

We may now substitute λ = 1 + δ into λ3 − 4λ2 + 5λ − 2 − 2 = 0 and this becomes

(1 + δ)3 − 4(1 + δ)2 + 5(1 + δ) − 2 = 2
=⇒ 1 + 3δ + 3δ 2 + δ 3 − 4 − 8δ − 4δ 2 + 5 + 5δ − 2 = 2
=⇒ δ 3 − δ 2 = 2 .

From (3.22), we get

δ 3 = c3 3 + 3c2 c2 4 + O
1
1

5

δ 2 = c2 2 + 2c1 c2 3 + 2c1 c3 + c2
1
2

Matching powers of

4

+O

in (3.23) gives the system of equations

− c2 = 1,
1
c3 − 2c1 c2 = 0,
1
3c2 c2 − 2c1 c3 − c2 = 0
1
2

which has a solution of

c1 = ±i,

c2 =

−1
,
2

53

c3 =

5i
8

5

.

(3.23)

and thus, the two complex roots of (3.21) are given by
(kc∆t)2 5(kc∆t)3 i
+ O ∆t4
2
8
3
±ikc∆t − 11(kc∆t) i + O ∆t4
=e
24
11(kc∆t)3 i
ω
=⇒ ei˜ ∆t = e±ikc∆t −
+ O ∆t4 .
24
λ = 1 ± ikc∆t −

Thus, we have that i˜ ∆t matches ±ikc∆t to the third order which implies that
ω

ω 2 = (kc)2 + O ∆t2
˜

and the phase error is of second order.

3.2.3

Fully Discrete Non-dissipative Scheme

Consider the non-dissipative scheme employing the midpoint quadrature rule for spatial
discretization. The only modification to equation (3.18) will be that the integral is replaced
ω
with an infinite sum. Let un = uei(kj∆x−n˜ ∆t) , where now ω denotes the fully discrete
ˆ
˜
j

temporal frequency, yielding
∞

ω
ω
ω
e−i˜ ∆t + ei˜ ∆t ueikj∆x−i˜ (n∆t) = α
ˆ

∆x/2

=−∞ −∆x/2

54

ω
e−α|x−y| ueiky−i˜ (n∆t) dy.
ˆ

ω
After dividing through by ueikj∆x−i˜ n∆t , this becomes
ˆ

ω
e−i˜ ∆t

ω
+ ei˜ ∆t

∞

αeik(m−j)∆x

=
m=−∞
∞

Ωm

αeik(m−j)∆x

=⇒ 2 cos (˜ ∆t) =
ω
m=−∞

e

−α|xj −y|

Ωm

e

dy

−α|xj −y|

dy.

The integration produces

−α|xj −y|

Ωm

e

dy =



2

 α 1 − e(−ν/2) ,

j=m


 2 −α|j−m|∆x
ν
 e
α
, j = m.
sinh
2

This then gives that



2 cos (˜ ∆t) = 2 1 − e(−ν/2) +
ω

eik(m−j)∆x e−α|j−m|∆x sinh
m=j

ν 
.
2

After dividing by two and centering the sum about m = j, this now becomes

cos (˜ ∆t) = 1 − e(−ν/2) + sinh
ω

ν
2

eikm∆x e−α|m|∆x
m=0


= 1 − e(−ν/2) + sinh
= 1 − e(−ν/2) + sinh
= 1 − e(−ν/2) + sinh

ν 
2
ν
2
ν
2

−∞

e(ik+α)∆xm +

m=−1
∞
e(−ik−α)∆xm

m=1
e(−ik−α)∆x

1 − e(−ik−α)∆x

55

+

∞
m=1
∞

+


e(ik−α)∆xm 
e(ik−α)∆xm

m=1
e(ik−α)∆x

1 − e(ik−α)∆x

.

Finding a common denominator then yields

ν
2

cos (˜ ∆t) = 1 − e(−ν/2) + sinh
ω
=

1 − e(−ν/2)

ν
+ sinh
2

e(−ik−α)∆x − 2e−2α∆x + e(ik−α)∆x
1 − e(−ik−α)∆x − e(ik−α)∆x + e−2α∆x
cos (k∆x) − e−ν
cosh (ν) − cos (k∆x)

.

Further algebraic manipulation shows that this becomes

cos (˜ ∆t) =
ω

1 − e(−ν/2)

ν
+ 2 sinh
2

cos (k∆x)eν − 1
1 + e2ν − 2 cos (k∆x)eν

=

1 − e(ν/2)

=

1 − e(ν/2)

ν
+ 2 sinh
2

cos (k∆x)eν − 1

ν
ν
+ 2 sinh
+ 2 sinh
2
2

= 1 − e(ν/2) + 2 sinh

1 + e2ν − 2 cos (k∆x)eν

1 + e2ν − 2 cos (k∆x)eν + cos (k∆x)eν − 1
1 + e2ν − 2 cos (k∆x)eν
1 − cos (k∆x)e−ν

ν
2

1 + e−2ν − 2 cos (k∆x)e−ν

.

As in the semi-discrete analysis, we study the phase for a fixed frequency k. Some
additional algebraic manipulation, including the double-angle formulas for cos and cosh,
gives

cos (˜ ∆t) = 1 − e(ν/2) + 2 sinh
ω

=1−

ν
2

1 − cos (k∆x)e−ν
1 − 2 cos (k∆x)e−ν + e−2ν

e(ν/2) 1 − 2 cos (k∆x)e−ν + e−2ν + 2 sinh (ν/2) 1 − 2 cos (k∆x)e−ν
1 − 2 cos (k∆x)e−ν + e−2ν
e(−ν/2) + e(−3ν/2) cos (k∆x) − e(−3ν/2) − e(−ν/2)

=1+

1 − 2 cos (k∆x)e−ν + e−2ν
e(ν/2) + e(−ν/2) cos (k∆x) − e(−ν/2) − e(ν/2)

=1+

eν − 2 cos (k∆x) + e−ν

56

.

This then becomes

cos (˜ ∆t) = 1 + 2 cosh
ω

ν
2

4 sinh2 (ν/2) + 4 sin2 (k∆x/2)

= 1 − 2 sin2

=⇒

cos (k∆x) − 1
2 cosh (ν) − 2 cos (k∆x)

= 1 + 2 cosh

=⇒ −

ν
2

k∆x
2

cos (˜ ∆t)
ω
−1
=
+ sin2
2
2

cos (k∆x) − 1

k∆x
2

cosh (ν/2)
2 sinh2 (ν/2) + 2 sin2 (k∆x/2)
cosh (ν/2)
2 sinh2 (ν/2) + 2 sin2 (k∆x/2)

1 cos (˜ ∆t)
ω
−
= sin2
2
2

k∆x
2

2 sinh2 (ν/2) + 2 sin2 (k∆x/2)

ω ∆t
˜
2

k∆x
2

2 sinh2 (ν/2) + 2 sin2 (k∆x/2)

=⇒ sin2

= sin2

cosh (ν/2)
cosh (ν/2)

.

Since we are interested in implicit methods, the relevant inequality is ∆x < c∆t, and so
consider the limit of the phase error for ∆x, ∆t → 0, while leaving ν as a small, fixed
quantity. The error will be better understood if put it in a form similar to equation (3.20).
Therefore, define

z=

kc∆t
,
2

s1 =

2
ν
sinh
,
ν
2

57

s2 =

2
sin
k∆x

k∆x
.
2

(3.24)

Then the fully discretized phase error can be given as

sin2

ω ∆t
˜
2

=

1 2
sin
2
sinh2

kνc∆t
√
2 2

ν
2
kνc∆t
√
2 2

cosh

ν
+ sin2
2

(zs2 )2
ν
2 s2 + 2(zs2 )2
1
ν
1
= z 2 cosh
˜
2 1 + 2˜2
z

= cosh

where

sin (k∆x/2)
1
s
= sin
z = 2z =
˜
s1
2 sinh (ν/2)
2

k∆x
ν
csch
.
2
2

(3.25)

Notice that since s1 = 1 + O ν 2 and s2 = 1 + O (k∆x)2 , that the semi-discrete limit
z ∼ z is recovered by letting ν and ∆x tend to zero such that ∆t is held fixed.
˜
The relative phase error will now depends on ν, ∆x and ∆t with

Φ(ν, ∆x, ∆t) =

ω
˜
−1 .
kc

Following the exact algebra and expansions as given in section 3.2.1, we have

ω = cosh (ν/2)kc 1 −
˜
= cosh

(kc∆t)2
4

+

cosh (ν/2)kc∆t
(kc∆t)2
1−
2
4
3∆t

ν
(kc∆t)2 (cosh (ν/2)(kc∆t))
kc 1 −
+
+ O ∆t4
2
4
24

58

.

3

+ O ∆t4

Using the expansion

cosh

ν
2

=1+

ν2
+ O ν4 ,
8

it becomes easy to see that the error is now O ∆t2 + ν 2 . If for a fixed wave number k, we
let ∆t and ∆x → 0 in such a way that ν is held constant, we see that the phase error is
non-zero in the static limit

˜
arcsin z
lim Φ(ν, ∆x, ∆t) =

∆x→0

cosh (ν/2)
1 + 2˜2
z
kc

−1 =

cosh (ν/2)
−1 ,
s1

which has been verified numerically for a range of ν. Additionally, letting ν → 0 produces
zero phase error, but this limit is not realizable, as it implies that either ∆x = 0 or ∆t = ∞.
The fully discrete phase error is shown for the second order method in Figure 3.1, for
several values of ν. Clearly the phase error is reduced as ν is reduced, which in the context
of the plot corresponds to spatial refinement.

3.2.4

Fully Discrete Dissipative Schemes

We now repeat the same ideas for the dissipative schemes. Begin once again with the ansatz
ω
un = uei(kj∆x−n˜ ∆t) , which when substituted into the first order scheme gives
ˆ
j

ω ˆ
ω
e−i˜ ∆t ueikj∆x−i˜ (n∆t)

∞

∆x/2

=α
=−∞ −∆x/2

59

ω
ω
e−α|x−y| ueiky−i˜ (n∆t) 2 − ei˜ ∆t dy.
ˆ

Relative phase error
0.9
0.8

c∆t=∆x/2
c∆t=∆x
c∆t=2∆x

0.7
0.6

Phi

0.5
0.4
0.3
0.2
0.1
0
0

0.1

0.2

0.3

0.4

0.5
k∆x/π

0.6

0.7

0.8

0.9

1

Figure 3.1: Relative Phase Error for the Non-dissipative Scheme
"For interpretation of the references to color in this and all other figures, the reader is referred
to the electronic version of this dissertation."

60

Repeating the integration techniques used in the previous section, we arrive at

ω
ω
e−i˜ ∆t = 2 − ei˜ ∆t σ

where

σ = 1 − e(ν/2) + 2 sinh

ν
1 − cos (k∆x)e−ν
.
2 1 − 2 cos (k∆x)e−ν + e−2ν

(3.26)

Using the exact algebra as was given above, we have in terms of the expressions (3.24) and
(3.25)

σ = 1 − cosh

sin2 (k∆x/2)
4˜2
z
ν
ν
= 1 − cosh
2 (ν/2) + sin2 (k∆x/2)
2 sinh
2 1 + 4˜2
z

ω
ω
=⇒ e−iˆ ∆t = 2 − eiˆ ∆t

1 − cosh

4˜2
z
ν
2 1 + 4˜2
z

.

(3.27)

ω
ω
Letting λ = eiˆ ∆t and then multiplying (3.27) by eiˆ ∆t yields the quadratic equation

σλ2 − 2σλ + 1 = 0.

This quadratic equation has solution λ = 1 ± i

1
σ

− 1. Here,

1
1 + 4˜2
z
−1=
−1
σ
1 + 4˜2 − cosh (ν/2)4˜2
z
z
=

4˜2 cosh (ν/2)
z
1 + 4˜2 − cosh (ν/2)4˜2
z
z

61

.

Therefore,
4˜2 cosh (ν/2)
z

1
−1=
σ

1 + 4˜2 − cosh (ν/2)4˜2
z
z

= 2˜
z

= 2˜
z

cosh (ν/2)
1 + 4˜2 − cosh (ν/2)4˜2
z
z
cosh (ν/2)
1 − 8˜2 sinh (ν/4)
z

.

This results in

ω
λ = eiˆ ∆t = 1 ± i sin

k∆x
ν
csch
2
2

cosh (ν/2)
1 − 8 sinh2 (ν/4)˜2
z

.

(3.28)

We now take the following expansions for each of the terms that make up the last term in
equation (3.28):

sin

k∆x
2

ν
2
ν
cosh
2

k∆x (k∆x)3
−
+ O ∆x5 ,
2
48

2
ν
+
+ O ν3 ,
ν 12
ν
= 1 − 1 + cosh
2
cosh (ν/2) − 1
ν2
=1+
+ O ν4 = 1 +
+ O ν4 ,
2
8
ν 2
ν 2
1 − 8 sinh2
z = 1 + 4 sinh2
˜
z + O ν2 .
˜
4
4
csch

and

=
=

62

Substituting all of these in to equation (3.28) gives

ω
λ = eiˆ ∆t =

k∆x
ν

1+

ν2
8

= 1 ± ikc∆t 1 +

1 + 4 sinh2

ν 2
z + O ∆x2 + ν 2
˜
4

ν2
ν
+ sinh2
(kc∆t)2 + O ∆x2 + ν 2
8
4

.

We then take the log of both sides and expand the right hand side using the Taylor
expansion for log (1 + x) to get
ν
ν2
+ sinh2
(kc∆t)2 + O ∆x2 + ν 2
iˆ ∆t = ±ikc∆t 1 +
ω
8
4

Dividing through by i∆t will give that the error is first order.

63

Chapter 4
Analysis of the Solution on a Bounded
Domain
We now repeat the analysis of these schemes on a bounded domain.

4.1
4.1.1

Consistency
Non-dissipative Scheme

First, consider the purely dispersive solver (2.6) with β = 2 on the region [−L, L] with either
zero Dirichlet or Neumann boundary conditions. The truncation error is now given as

τn =

1
(c∆t)2

v n+1 + v n−1 − J [v n ] ,

64

where
L

J [v n ] = α

e−α|x−y| v n (y) dy + c1 (x)
˜

−L
L

+ c2 (x)
˜

L

e−α(L−y) v n (y) dy

(4.1)

−L

e−α(L+y) v n (y) dy .

−L

Zero Dirichlet boundary conditions will give

c1 (x) = −
˜

sinh(α(L + x))
,
sinh(2αL)

c2 (x) = −
˜

sinh(α(L − x))
sinh(2αL)

while zero Neumann boundary conditions give

c1 (x) =
˜

cosh(α(L + x))
,
sinh(2αL)

c2 (x) =
˜

cosh(α(L − x))
.
sinh(2αL)

This can be written as the integral of a single function Q(x, y) multiplied by v n with

Q(x, y) = α e−α|x−y| + c1 (x)e−α(L−y) + c2 (x)e−α(L+y) .
˜
˜

and due to the negative exponentials, it is easy to see that

|Q(x, y)| < 2,

∀(x, y) ∈ [−L, L].

Observe that a direct consequence of the homogeneous Dirichlet boundary conditions is
that at x = ±L, v n and all of its even spatial derivatives will vanish. That is,

lim (∂x )2m v(x, t) = lim

x→±L

1

x→±L c2m

(∂t )2m v(x, t) =

65

1
c2m

(∂t )2m v(±L, t) = 0.

The order of the limits and the derivatives can be swapped as long as the 2m-th derivatives
of v n exist and are continuous.
Similarly, for homogenous Neumann boundary conditions, the odd derivatives will disappear as

lim (∂x )2m+1 v(x, t) = lim

1

x→±L c2m

x→±L

(∂t )2m vx (x, t) =

1
c2m

(∂t )2m vx (±L, t) = 0.

Now, we perform integration by parts on J [v n ] = J+ [v n ] + J− [v n ] + c1 (x)K+ [v n ] +
˜
c2 (x)K− [v n ], where
˜
L x

J± = α

0

e−αz v n (x ± z) dz
L x

= −v n (x ± z)e−αz 0

L x

+

0

n
(±1)vx (x ± z)e−αz dz

.
.
.


k

= −
j=0
L x

−

0

±1
∂x
α
±1
∂x
α

j

L x
v n (x ± z)e−αz 
0
k+1

v n (x ± z)e−αz dz.

66

Similarly, the boundary terms are
L

K− = α

e−α(L+y) v n (y) dy

−L
L

L

= −v n (y)e−α(L+y)

−L

+
−L

n
vx (y)e−α(L+y) dy

.
.
.


k

= −
j=0
L

+
−L

j

1
∂x
α

L
v n (y)e−α(L+y) 
−L

1
∂x
α

k+1

v n (y)e−α(L+y) dy

and
L

K+ = α

e−α(L−y) v n (y) dy

−L

= v n (y)e−α(L−y)

L

L

−

−L

−L

n
vx (y)e−α(L−y) dy

.
.
.


k

=
j=0
L

+
−L

−1
∂x
α

j

L
v n (y)e−α(L−y) 
−L

−1
∂x
α

k+1

v n (y)e−α(L−y) dy.

Taking this expansion out to k = 5 and after grouping by odd and even derivatives at

67

the boundary points, we arrive at

1 n
1 n
vxx (−L) +
vxxxx (−L)
α2
α4
1 n
1 n
1 n
vx (−L) +
v (−L) +
vxxxxx (−L)
+ δ+
3 xxx
α
α
α5
1 n
1 n
vxx (L) +
vxxxx (L)
+ γ− v n (L) +
α2
α4
1 n
1 n
1 n
vx (L) +
vxxx (L) +
vxxxxx (L)
+ γ+
α
α3
α5

J [v n ] = δ− v n (−L) +

+ 2 v n (x) +
+

L

1
α6

0

1 n
1 n
v (x) +
vxxxx (x)
2 xx
α
α4

n
Q(x, y)vxxxxxx (y) dy

with

δ± = c2 (x) ± e−2αL c1 (x) ± e−α(L+x) ,
˜
˜
γ± =

c1 (x) − e−2αL c2 (x) − e−α(L−x) .
˜
˜

For zero Dirichlet boundary conditions,

sinh(2αL)δ+ = −sinh(α(L − x)) − e−2αL sinh(α(L + x)) + sinh(2αL)e−α(L+x)
=

−eα(L−x) + eα(−L+x) − eα(−L+x) + eα(−3L−x) + eα(L−x) − eα(−3L−x)
2

=0

68

while

sinh(2αL)δ− = −sinh(α(L − x)) + e−2αL sinh(α(L + x)) − sinh(2αL)e−α(L+x)
=

−eα(L−x) + eα(−L+x) + eα(−L+x) − eα(−3L−x) − eα(L−x) + eα(−3L−x)
2

= eα(−L+x) − eα(L−x) = 2sinh(α(x − L)).

Additionally,

sinh(2αL)γ+ = sinh(α(L + x)) + e−2αL sinh(α(L − x)) − sinh(2αL)e−α(L−x)
=

eα(L+x) − eα(−L−x) + eα(−L−x) − eα(−3L+x) − eα(L+x) + eα(−3L+x)
2

=0

and

sinh(2αL)γ− = −sinh(α(L + x)) + e−2αL sinh(α(L − x)) − sinh(2αL)e−α(L−x)
=

−eα(L+x) + eα(−L−x) + eα(−L−x) − eα(−3L+x) − eα(L+x) + eα(−3L+x)
2

= −eα(L+x) + eα(−L−x) = −2sinh(α(L + x)).

However, it was already shown that all even derivatives of v n are zero at the boundary.
Thus, all of the boundary terms vanish identically, and

J [v n ] = 2v n +

L
2 n
2 n
1
n
vxx +
vxxxx +
Q(x, y)vxxxxxx (y) dy.
2
4
6 −L
α
α
α

(4.2)

A comparison of (4.2) and (3.2) shows that J [v n ] and I[v n ] differ by O ∆t6 . Thus, for

69

all cases on a bounded domain with Dirichlet boundary conditions, J [v n ] may be substituted
for I [v n ] into the free space consistency proof for the respective solver without having an
effect on the order of the method.
For zero Neumann boundary conditions,

sinh(2αL)δ+ = cosh(α(L − x)) + e−2αL cosh(α(L + x)) + sinh(2αL)e−α(L+x)
eα(L−x) + eα(−L+x) + eα(−L+x) + eα(−3L−x) + eα(L−x) − eα(−3L−x)
=
2
= eα(L−x) + eα(−L+x) = 2cosh(α(L − x))

while

sinh(2αL)δ− = cosh(α(L − x)) − e−2αL cosh(α(L + x)) − sinh(2αL)e−α(L+x)
=

eα(L−x) + eα(−L+x) − eα(−L+x) − eα(−3L−x) − eα(L−x) + eα(−3L−x)
2

= 0.

Additionally,

sinh(2αL)γ+ = −cosh(α(L + x)) − e−2αL cosh(α(L − x)) − sinh(2αL)e−α(L−x)
=

−eα(L+x) − eα(−L−x) − eα(−L−x) − eα(−3L+x) − eα(L+x) + eα(−3L+x)
2

= −eα(−L−x) − eα(L+x) = −2cosh(α(L + x))

70

and

sinh(2αL)γ− = cosh(α(L + x)) − e−2αL cosh(α(L − x)) − sinh(2αL)e−α(L−x)
=

eα(L+x) + eα(−L−x) − eα(−L−x) − eα(−3L+x) − eα(L+x) + eα(−3L+x)
2

= 0.

However, it was already shown that all odd derivatives of v n are zero at the boundary.
Thus, all of the boundary terms vanish identically, and

J [v n ]

=

2v n

L
2 n
1
2 n
n
v +
v
+
Q(x, y)vxxxxxx (y) dy.
+
2 xx
4 xxxx
6 −L
α
α
α

Once again, we may freely substitute J [v n ] for I [v n ] in any of the freespace consistency
proofs without changing the order.

4.1.2

Dissipative Schemes

For the dissipative schemes (2.3) and (2.4), we can show consistency based on results obtained
in previous sections. From section 3.1.1.2, we have that the first order scheme is given by

1
un+1 (x) = I [un ] − I un−1
2

(4.3)

and the second order scheme is given by

5
1
un+1 (x) = I [un ] − I un−1 + I un−2 .
4
4

71

(4.4)

We can adapt these situations to the case of a finite mesh as

1
un+1 (x) = J [un ] − J un−1
2

(4.5)

5
1
un+1 (x) = J [un ] − J un−1 + J un−2
4
4

(4.6)

and

respectively, where J is as defined in (4.1). From section 3.1.1.2, we have that equations
(4.3) and (4.4) are first and second order, respectively. From section 4.1.1, we have that
J [u] = I [u] + O ∆t6 . Thus, (4.5) and (4.6) must similarly be of first and second order.

4.2

Stability

In this section, the focus is on the incorporation of the boundary matrix, such that At = Ap +
Ah in equation (3.5). We will treat the matrix Ah , given by equation (3.8), perturbatively
and look for a bound on the eigenvalues of At , based on λ + δ, where λ is an eigenvalue of
Ap and δ is an eigenvalue of Ah . Stability will follow if it can be shown that |µ| ≤ 1, where
µ is now given by

µ2 − (λ + δ)µ +

λ+δ
= 0 and
2

5
1
µ3 − (λ + δ)µ2 + (λ + δ)µ − (λ + δ)µ = 0
4
4

72

in the respective dissipative schemes and

µ2 − (λ + δ)µ + 1 = 0

in the non-dissipative scheme. While stability will hold for the non-dissipative scheme if
λ + δ ∈ [−2, 2], the interval is [−2/3, 2] for the first order dissipative scheme and [−2/5, 2]
for the second order dissipative scheme. The region for the second order dissipative scheme
was found numerically while the other regions may be found by solving a quadratic equation.
Lemma 6. λ + δ ≤ 2.
Proof. By the Gerschgorin theorem, the maximal eigenvalue of Ap + Ah is bounded above
by the maximal row sum of Ap + Ah . For zero Dirichlet boundary conditions,
L

(Ap )ij =
j

αe−α|xi −y| dy = 2 − e−α(xi +L) − e−α(L−xi ) and

−L

sinh(α(L + xi )) L −α|L−y|
(Ah )ij = −
αe
dy
sinh(2αL)
−L
j
sinh(α(L − xi )) L −α|−L−y|
αe
dy
sinh(2αL)
−L
−sinh(α(L + xi )) − sinh(α(L − xi ))
= 1 − e−2αL
.
sinh(2αL)

−

Thus,

(Ap + Ah )ij = 2 − e−α(xi +L) − e−α(L−xi )
j

+ 1 − e−2αL

−sinh(α(L + xi )) − sinh(α(L − xi ))
< 2.
sinh(2αL)

73

In a similar manner, zero Neumann boundary conditions will give
L

(Ap )ij =
j

(Ah )ij =
j

αe−α|xi −y| dy = 2 − e−α(xi +L) − e−α(L−xi ) and

−L

cosh(α(L + xi )) L −α|L−y|
αe
dy
sinh(2αL)
−L

cosh(α(L − xi )) L −α|−L−y|
αe
dy
sinh(2αL)
−L
cosh(α(L + xi )) + cosh(α(L − xi ))
= 1 − e−2αL
sinh(2αL)
+

=

e−α(L+xi ) + eα(L+xi ) + eα(xi −L) + eα(L−xi )
2sinh(2αL)
−

e−α(3L+xi ) − eα(xi −L) − eα(xi −3L) − e−α(L+xi )
2sinh(2αL)

e−αL e2αL − e−2αL
=

eαxi + e−αxi

2sinh(2αL)

= 2e−αL cosh(αxi ).

Thus,

(Ap + Ah )ij = 2 − e−α(xi +L) − e−α(L−xi ) + 2e−αL cosh(αxi )
j

= 2 − e−α(xi +L) − e−α(L−xi ) + e−α(xi +L) + e−α(L−xi ) = 2.

Since we already know that λ > 0, we must ensure that −2/3 < δ < 0 in the first order
scheme and −2/5 < δ < 0 in the second order dissipative scheme to ensure stability. Below,
we will make use of the Bauer-Fike theorem.

74

Theorem 7 (Bauer-Fike Theorem [11]). Let B be a diagonalizable matrix, and Q be the
non-singular eigenvector matrix such that A = QΛQ−1 , where Λ is the diagonal matrix
containing the eigenvalues of B. Then for the matrix B + E, with eigenvalues λ + δ, and
where E is an arbitrary perturbation matrix,

|δ| ≤ ||Q||2 ||Q−1 ||2 ||E||2 .

We will apply this theorem with B = Ap , which is a real symmetric positive definite
matrix, and therefore we have that ||Q||2 ||Q−1 ||2 = 1 since Q is orthogonal. Also, set
E = Ah and thus

|δ| ≤ ||Ah ||2 .

Therefore, attention will be restricted to the norm of Ah .
Notice that if we define the vectors

vd = (sinh (α(L − x1 )), . . . sinh (α(L − xN )))T ,
vn = (cosh (α(L − x1 )), . . . cosh (α(L − xN )))T ,
h = e−α(L+x1 ) , . . . e−α(L+xN )

then the matrix Ah (3.8) can be expressed as

Ah = −2a0

sinh (ν/2)
(v · h + (Jvi ) · (Jh)),
sinh (2αL) i

i ∈ {d, n}

(4.7)

where a0 is given by (3.9) and the matrix J is the counter-identity matrix, containing ones
75

along the anti-diagonal. It has the effect of reversing the elements of a vector that it acts
upon, and for an equally spaced symmetric mesh is equivalent to replacing xj with −xj .
Since v · h and (Jv) · (Jh) are both rank one matrices, Ah is a matrix of at most rank two.
Theorem 8. For the centrosymmetric rank-two matrix with n distinct entries which are
negative and of decreasing magnitude less than one, the separation between the eigenvalues
is bounded by 2an .
The proof of this theorem will make use of the following lemmas.
Lemma 9. For these centrosymmetric matrices, Ah (1, j)2 + Ah (1, n − j + 1)2 = Ah (1, 1) ∗
Ah (j, j) + Ah (1, n) ∗ Ah (j, n − j + 1).
Proof. Letting ν =

dx
sinh((n − i + 0.5)ν) −(j−0.5)ν sinh((i − 0.5)ν) −(n−j+0.5)ν
, Ah (i, j) =
e
+
e
.
c∆t
sinh(nν)
sinh(nν)

This gives then that

sinh2 (nν)(Ah (1, j)2 + Ah (1, n − j + 1)2 ) =
sinh(−0.5ν)e(j−n−0.5)ν + sinh((0.5 − n)ν)e(0.5−j)ν

2

sinh(−0.5ν)e(0.5−j)ν + sinh((0.5 − n)ν)e(j−n−0.5)ν

2

+

= 4 + e(2j−2)ν + e(2n−2j)ν − 4e−ν + e−2jν + e(2−2j)ν + e(2j−2n)ν − 4e(1−2n)ν
+ 4e−2nν + e(2j−4n)ν − 4e(1−2j)ν − 4e(2j−2n−1)ν + e(2−2n−2j)ν + e(2j−2n−2)ν

76

which is equal to

sinh(−0.5ν)e(0.5−n)ν + sinh((0.5 − n)ν)e−0.5ν
· sinh((0.5 − j)ν)e(j−n−0.5)ν + sinh((j − n − 0.5)ν)e(0.5−j)ν
+ sinh(−0.5ν)e−0.5ν + sinh((0.5 − n)ν)e(0.5−n)ν
· sinh((0.5 − j)ν)e(0.5−j)ν + sinh((j − n − 0.5)ν)e(j−n−0.5)ν
= sinh2 (nν)(Ah (1, 1) ∗ Ah (j, j) + Ah (1, n) ∗ Ah (j, n − j + 1)).

For the case of zero Neumann boundary conditions,

A(i, j) =

cosh((n − i + 0.5)ν)
cosh((i − 0.5)ν)
exp((j − 0.5)ν) +
exp((n − j + 0.5)ν),
sinh(nν)
sinh(nν)

which will consistently change signs in the algebra above and the result will be the same.
The corollary seen from this is that

n
2
i=1 ai

a Trace(An ) + an SkewTrace(An )
.
= 1
2

Lemma 10. For the given matrix, SkewTrace(An )2 = Trace(An )2 −4Mn , where Mn consists
of the principal 2-minors of An .
Proof. Proof of theorem by induction:


a1 a2 
Base Case: Consider the two by two matrix A2 = 
.
a2 a1

Trace(A2 )2 − 4Mn = (2a1 )2 − 4 a2 − a2
1
2
= 4a2 = SkewTrace(A2 )2 .
2

77

Induction Step:
Assume this true for any such matrix of size n. Extend the matrix to An+2 of size n + 2 by
adding on an additional first and last row and column. Due to the nature of the extension
from An to An+2 , Trace(An+2 ) = Trace(An ) + 2An+2 (1, 1). Thus,

Trace(An+2 )2 = Trace(An )2 + 4An+2 (1, 1)Trace(An ) + 4An+2 (1, 1)2 .

(4.8)

Mn can be extended to Mn+2 by adding in the principal 2-row minors that have an entry
in the first or last row. Thus,

Mn+2 = Mn + An+2 (1, 1)2 − An+2 (1, n + 2)2
n

+ 4An+2 (1, 1)Trace(An ) − 2

(4.9)

An+2 (1, i)2 .

i=2

Thus,

Trace(An+2 )2 − 4Mn+2
= Trace(An )2 + 4An+2 (1, 1)2 + 4An+2 (1, 1)Trace(An )
− 4 Mn + An+2

(1, 1)2

− An+2

(1, n + 2)2

n−1

+ 2An+2 (1, 1)Trace(An ) − 2

An+2 (1, i)2

i=2

= Trace(An )2 − 4Mn + 4An+2 (1, n + 2)2 − 4An+2 (1, 1)Trace(An ) + 8

n
i=2

78

An+2 (1, i)2

From the base case, this is now

SkewTrace(An

)2

+ 4An+2

(1, n + 2)2

n

− 4An+2 (1, 1)Trace(An ) + 8

An+2 (1, i)2

i=2

= SkewTrace(An )2 + 4An+2 (1, n + 2)2 + 4An+2 (1, n + 2)SkewTrace(An )
= (SkewTrace(An ) + 2An+2 (1, n + 2))2
= SkewTrace(An+2 )2 .

Proof. Proof of theorem 8 by induction:
Base case 1:
Consider the case of the one by one centrosymmetric matrix [a]. This obviously has one
eigenvalue a and there is no separation to consider.
Base case 2:


a1 a2 
Consider the case of the two by two centrosymmetric rank-two Hankel matrix 

a2 a1
which obviously has eigenvalues λ = a1 ± a2 , which have a difference of 2a2 = 2an .

Induction Step:
Assume this is true for any matrix An of size n. Extend the matrix to An+2 of size n + 2
by adding on an additional first and last row and column. The characteristic polynomial is
given as
2n

ci λi .

P (λ) =
i=0

79

(4.10)

Since this is the characteristic polynomial, c2n = 1. Further, because the matrix is of ranktwo, ci = 0 for i < 2n − 2. Also, c2n−1 = Trace(An+2 ) and c2n−2 is the sum of the principal
2-row minors. Solving the resultant quadratic equation yields that the eigenvalues of the
matrix are

1
2

Trace(An+2 ) +

Trace(An+2 )2 − 4Mn+2

(4.11)

where Mn+2 denotes the sum of the principal 2-row minors. Since An and An+2 are real and
symmetric, all the eigenvalues must be real, which implies that Trace(An+2 )2 ≥ 4Mn+2 .
By lemma 10, the eigenvalues for An+2 are Trace(An+2 ) ± SkewTrace(An+2 ).

From the above theorem, the eigenvalues of Ah are given by

Trace(Ah ) ± Trace(JAh )
.
2

Furthermore, it is clear from equation (3.8) that Aij = Aji , giving that Ah is symmetric.
Therefore, the 2-norm of the matrix is equal to its spectral radius

Trace(Ah ) + Trace(JAh )
2
sinh (ν/2)
(v · h + (Jv) · (Jh) + (Jv) · h + v · (Jh)) .
=
sinh (2αL)

||Ah ||2 = −

These dot products are each given by a geometric series determined by the boundary

80

conditions. For Dirichlet boundary conditions,
N

v · h = (Jv) · (Jh) =

−α(L+xj )

e

sinh α(L − xj )

j=1

=

=

1
2

N

e

−2αxj

− e−2αL

j=1

1
2

sinh (N ν)
− N e−N ν
sinh (ν)

and
N

(Jv) · h = v · (Jh) =

e

−α(L+xj )

sinh α(L + xj )

j=1

=

=

1
2
1
2

N

−2α(L+xj )

1−e
j=1

N − e−N ν

sinh (N ν)
sinh (ν)

,

where we recall that 2αL = N ν. Thus,

|δ| ≤ ||Ah || = a0 1 − e−N ν
= a0 1 − e−N ν

sinh (ν/2)
sinh (N ν)

sinh (N ν)
+N
sinh (ν)

sinh (ν/2)
sinh (ν/2)
+N
sinh (ν)
sinh (N ν)

81

.

Neumann boundary conditions give
N

v · h = (Jv) · (Jh) =

e

−α(L+xj )

cosh α(L − xj )

j=1

=

=

1
2

N

−2αxj

e

+ e−2αL

j=1

1
2

sinh (N ν)
+ N e−N ν
sinh (ν)

and
N

(Jv) · h = v · (Jh) =

e

−α(L+xj )

cosh α(L + xj )

j=1

=

=

1
2
1
2

N

−2α(L+xj )

1+e
j=1

N + e−N ν

sinh (N ν)
sinh (ν)

.

Thus,

|δ| ≤ ||Ah || = a0 1 + e−N ν
= a0 1 + e−N ν

sinh (ν/2)
sinh (N ν)

sinh (N ν)
+N
sinh (ν)

sinh (ν/2)
sinh (ν/2)
+N
sinh (ν)
sinh (N ν)

.

Figures 4.1 and 4.2 give numerical plots of these eigenvalue bounds versus ν = ∆x/(c∆t)
and N . The figure encompasses 1 ≤ N ≤ 75 and 0.01 ≤ ν ≤ 2. The white strips indicate
a region where the bound δ > −2/3 is not satisfied, and therefore indicates a potential for
instabilities to arise. The blue area indicates a region where −2/3 < δ < −2/5 indicating
stability of the first order scheme but not of the second order scheme. However, no such
82

instabilities have been observed in my numerical experiments thus far, and it is likely that
if a tighter bound were taken, this region of instability would not be present. It should
be emphasized that this stability region only applies to the dissipative schemes and has no
bearing on the second order purely dispersive scheme.

4.3

Numerical Results

In this section we present preliminary results for the convergence of the outlined numerical
schemes for a test case. The wave equation is computed on the unit interval, with wave
speed c = 2, up to a time t = T = 1, and zero Dirichlet conditions are imposed. For the
initial condition, we prescribe a Gaussian function
12(x − 0.5) 2
−
L
f (x) = e
,

g(x) = 0

which will be initially of negligible amplitude at the boundaries, and will have the exact
solution

−

v(x, t) =

e

12(x − 0.5 − ct) 2
12(x − 0.5 + ct) 2
−
L
L
+e
.
2

The discretization is made so that the ratio c∆t/∆x is held fixed at 10 in the first order
method, and 20 in the second order methods. Note that in both cases, this quantity is well
beyond the CFL stability limit imposed on standard finite difference schemes, and so this
method is an order of magnitude faster. The spatial step ∆x = 1/N is varied in the range
83

Lower Eigenvalue Bound for Dirichlet Boundary Conditions
−0.4

Number of Mesh Points

20
−0.45

40
60

−0.5

80
100

−0.55

120
140

−0.6

160
180
200

−0.65
50

100

150
200
100*! x/! t

250

300

Figure 4.1: Contour Plot of the Minimum Eigenvalue of the Dissipative Schemes for Zero
Dirichlet Boundary Conditions with Midpoint Quadrature

Lower Eigenvalue Bound for Neumann Boundary Conditions
−0.4

Number of Mesh Points

20
−0.45

40
60

−0.5

80
100

−0.55

120
140

−0.6

160
180
200

−0.65
50

100

150
200
100*! x/! t

250

300

Figure 4.2: Contour Plot of the Minimum Eigenvalue of the Dissipative Schemes for Zero
Neumann Boundary Conditions with Midpoint Quadrature

84

corresponding to N = 50 ∗ 2k , k = {1, 2, . . . , 8}, and ∆t is scaled accordingly. The error is
computed as

||v(x, t) − u(x, t)||L2 [0,1] .

The error for the first order dissipative scheme is shown in Figure 4.3. The convergence
is as expected, with the exception of the gradual decrease for larger values of ∆t. Similar
trends are observed in Figures 4.4 and 4.5, the second order schemes that are dissipative
and non-dissipative, respectively. Since our interest is in taking time steps larger than the
CFL condition, it is reassuring to see the error drop with increasing time step over some
nontrivial interval.

Convergence Study for the First Order Dissipative Scheme

1

10

Midpoint
Trapezoidal
!t1

0

||u−v||2

10

−1

10

−2

10 −4
10

−3

10

−2

!t

10

−1

10

Figure 4.3: Convergence of the First Order Dissipative Scheme with ν = ∆x/(c∆t) = 0.1.

85

Convergence Study for the Second Order Dissipative Scheme
2
10
Midpoint
Trapezoidal
1
!t2
10

||u−v||2

0

10

−1

10

−2

10

−3

10 −4
10

−3

10

−2

!t

10

−1

10

√
Figure 4.4: Convergence of the Second Order Dissipative Scheme with ν/ 2 = ∆x/(c∆t) =
0.05.

86

Convergence Study for the Second Order Non−Dissipative Scheme
2
10
Midpoint
Trapezoidal
1
10
!t2
0

||u−v||2

10

−1

10

−2

10

−3

10

−4

10 −4
10

−3

10

−2

!t

10

−1

10

√
Figure 4.5: Convergence of the Second Order Non-dissipative Scheme with ν/ 2 =
∆x/(c∆t) = 0.05.

87

4.4

Fast Implementation

All of the structure underlying both the matrices Ap and Ah creates opportunities for either
fast implementation or minimized storage. The matrix Ap (3.10), which corresponds to the
free space solution, can be written as

Ap = ˜0 I + ˜1 (L + U ),
b
b

where I is the standard identity matrix and L and U are strictly lower and upper triangular
matrices respectively. The entries of L and U have been defined as



0,


j≥i
Lij = Uji =

 −|i−j|ν
e

, j < i.
Consider the matrix vector multiplication Ap · x with x = [x1 , x2 , . . . , xN ]T . Then
¯
¯

Ap · x = ˜1 L¯ + ˜1 U x + ˜0 I x.
¯ b x b ¯ b ¯

a1 I x requires only that a1 be stored and multiplied by each entry in x, requiring precisely
¯
¯
N operations.
Both L¯ and U x can be quickly solved in a manner akin to Horner’s method. The
x
¯
structure of L will yield that

L¯ = [0, e−ν x1 , e−ν x2 + e−2ν x1 , e−ν x3 + e−2ν x2 + e−3ν x1 , . . . ].
x

88

Starting with the second entry, each entry can be found recursively based solely on the
previous entry using

L(x1 ) = 0,

L(xk ) = e−ν ∗ (L (xk−1 ) + xk−1 ) ,

k = {2, . . . , N }.

Similarly,

U (xN ) = 0,

U (xN −k ) = e−ν ∗ (U (xN −k+1 ) + xN −k+1 ) ,

k = {2, . . . , N }.

Each of these matrix multiplications will require N multiplications and N additions. Therefore, the matrix vector multiplication Ap · x requires storing only the quantities a0 , a1 , and
¯
e−ν as well as O(N ) computations.
The matrix Ah , corresponding to the boundary correction can also be stored as two O(N )
vectors as can easily be seen by equation (4.7). Additionally, since Ah may be defined as
a constant multiplied by the summand of two outer products, v · h and Jv · Jh, it may be
implemented fast as well. Note that multiplying the outer product by another vector x gives,
¯
for example,

(v · h)¯ = v · h · x = v · (h · x) = (h · x)v.
x
¯
¯
¯

Here, h · x is simply a dot product which requires only 2N − 1 operations and that scalar is
¯
simply multiplied by v. Thus, the matrix vector multiplication for Ah will require this to be
done twice, yielding a total of O(N ) operations.

89

Chapter 5
Non-Reflecting Boundary Conditions
A common problem that must be resolved when solving the wave equation is choosing how to
truncate the domain for numerical computation. For a wave propagating on a finite domain,
boundary conditions should be applied. While Dirichlet or Neumann conditions may be
used, a more practical approach may be to create an artificial non-reflecting boundary which
will allow the wave to flow out of the domain without re-entering [37]. One example of
such an approach is to use a radiation condition at infinity, which states that the waves
are outgoing. For Rd , d ∈ {1, 2, 3}, the most well-known such condition is the Sommerfeld
radiation condition [64], given by

lim r(d−1)/2 (ur − iku) = 0,

r→∞

(5.1)

where r is the radial coordinate, u is the scattered field and k is the wavenumber. For the
case of one dimensional time harmonic waves, u(x, t) = u(x)e−iωt , the Sommerfeld radiation
˜

90

condition (5.1) will exactly reduce to

ut ± cux = 0,

(5.2)

where the sign depends on the direction the wave is traveling. Both the accuracy [35] and the
stability [43] of solvers using various finite difference approximations implementing (5.2) have
been studied. In one dimension still, applying either of (5.2) is analogous to implementing
a one-way wave, that is, a wave that is allowed to propagate in only one direction [44, 68].
While a viable option in one dimension, radiation conditions may experience problems
in higher dimensions, however. The implementation of such conditions in two or three
dimensions has typically suffered from significant spurious reflections at the boundary [37].
Several attempts have been made to correct this, most notably [32, 33].
An alternative approach introduced by Berenger [15] is the idea of Perfectly Matched
Layers (PML), specifically designed for the context of solving Maxwell’s equation in two dimensions, wherein the computational domain is surrounded by a layer of absorbing material.
Berenger split the magnetic field Hz into directional components Hzx and Hzy and then
showed that if the electric conductivity σ and the magnetic conductivity σ∗ satisfied

σ
0

=

σ∗
,
µ0

then the interface between the computational domain and the absorbing material would be
perfectly matched allowing the wave to pass through into the absorbing medium without
reflection. The wave would then decay exponentially in the absorbing layer. Berenger would
later go on to adapt this method for the three dimensional case of Maxwell’s equations [16].

91

This idea of a perfectly matched layer has since been extended to wave problems in other
fields such as the linearized Euler equations [47], advective acoustics [1], shallow water waves
[58], etc.
The PML method proved numerically advantageous over other such methods for truncating domains, such as radiating boundary conditions and absorbing boundary conditions.
However, Abarbanel and Gottlied [3] showed that Berenger’s initial PML implementation
with field splitting was only weakly well-posed and thus when slightly perturbed, this numerical solution was subject to becoming divergent due to explosive modes. This splitting was
later shown to be unnecessary by Zhao and Cangellaris [72]. Many other studies, including
[2, 4, 7, 13, 12, 14, 30], have investigated the key properties of well-posedness and stability
of PML in its various schemes.
Another key tenet of PML that has been frequently studied is the method’s efficiency.
These implementations include further component splitting, coordinate stretching, the use
of absorbing boundaries in the absorbing medium [17], reductions in number of auxiliary
variables [42] and optimization of both the absorbing function and the absorbing material
grid [9, 18, 26]. While significant time has been invested into finding an optimal absorbing
function, the problem is still regarded as open.

5.1

Perfectly Matched Layers

Let v(x, t) be a true solution to the 1D wave equation (1.1). We wish to compute a numerical solution u(x, t) to this problem on the truncated computational domain [−L, L] by
introducing perfectly matched layers of length d. Boundary conditions at the end of the

92

perfectly matched layer may be enforced, but are unnecessary. Thus, we wish to compute

utt = c2 uxx ,
˜˜

t > 0, −L − d < x < L + d
˜

u(x, 0) = f (x),

(5.3)

ut (x, 0) = g(x)

u(±(L + d), t) = 0.

Here, we have chosen to enforce zero Dirichlet conditions at the end of the PML region for
convenience.
The new variable x will be defined using an analytic continuation of the wave equation
˜
and the Laplace transform. Let
∞

u(x, s) =
ˆ

0

e−st u(x, t)dt

be the Laplace transform of the solution to equation (1.1). Now, define x by the mapping
˜

d˜
x
σ(x)
=1+
,
dx
s

where σ(x) is some spatially varying damping rate, which is non-decreasing into the PML,
and nonzero only in the PML regions. Next, since

∂u
dx ∂u
s
∂u
=
=
,
∂x
˜
d˜ ∂x
x
s + σ(x) ∂x

93

we can take the Laplace transform of equation (5.3) (with f (x) = g(x) = 0), and get

s
∂x
s+σ

s 2 u = c2
ˆ
=⇒ s2 u =
ˆ

(sc)2
(s + σ)2

s
ux
ˆ
s+σ

uxx −
ˆ

σ
ux
ˆ
s+σ

s+σ 2
u = uxx − v
ˆ ˆ
ˆ
c

=⇒

where

v=
ˆ

σ
ux
ˆ
s+σ

=⇒

(s + σ)ˆ = σ ux .
v
ˆ

Now, using the inverse Laplace transform on the expressions for u and v gives the following
system of equations:

1
c2

utt + 2σut + σ 2 u = uxx − v,
vt + σv = σ ux .

(5.4)
(5.5)

Alternatively, we can write a second order integro-differential equation involving only u,
by solving for v using an integrating factor. Elimination then gives
t

1

utt + 2σut + σ 2 u = uxx − σ
2
c

94

0

ux (x, s)e−σ(t−s) ds.

(5.6)

5.1.1

Constant Damping Function

As a first attempt to understand this approach, consider the case where σ(x) is a piecewise
positive constant, with the wave equation satisfied for x < L, and the PML region for x > L;
then

σ(x) =





0,

x<L
.



σ > 0, x > L
 0
Then, σ(x) = σ0 H(x − L), where H(x) is the Heaviside step function, and therefore

σ (x) = σ0 δ(x − L).

5.1.1.1

(5.7)

Implementation

Making use of (5.7) in equation (5.6) produces

1
c2

utt + 2σut + σ 2 u = uxx − σ0 δ(x − L)

t
0

ux (x, s)e−σ(t−s) ds.

The delta function will not prove problematic for my method of solution since we intend to
integrate over space with the Green’s function.
We now discretize equation (5.4) in time, using the standard first order discretization to
get

1
(c∆t)2

un+1 − 2un + un−1 + 2σ∆t un+1 − un + (σ∆t)2 un+1 = un+1 − v n+1
xx

95

or

1
(c∆t)2

2un+1 − 5un + 4un−1 − un−2

+ 2σ∆t

3 n+1
1
u
− 2un + un−1
2
2

+ (σ∆t)2 un+1

= un+1 − v n+1
xx

by using the standard second order discretization.
The term v n+1 can be discretized in a number of ways. If we apply the integrating factor
to equation (5.5) and integrate over [tn , tn+1 ] using the trapezoidal rule, we get

v n+1 =

σ ∆t n+1
e−σ∆t v n + δ(x − L) 0
ux + e−σ∆t un
x
2

Now, we can apply the Helmholtz operator, which gives

L un+1 (x) = f un+1 , un , un−1 , v n+1

where

L[u] :=

∂2
∂x2

− α2 u

and α is dependent upon the discretization used.

96

.

The right hand side, depending on the desired order, is then given either by

f un+1 , un , un−1 , v n+1 = −α2 2un − un−1 − 2zun+1 + 2zun − z 2 un+1 + v n+1 or
f un+1 , un , un−1 , un−2 , v n+1 =
−

α2
5un − 4un−1 + un−2 − 3zun+1 + 4zun − zun−1 − z 2 un+1 + v n+1
2

with

z = σ∆t.

The Green’s function is given by

G(x|y) =

−1 −α|x−y|
e
2α

and upon inverting the Helmholtz operator, we have either

un+1 =α

∞

e−α|x−y| un (y) −

−∞
∞

−α
L

un−1 (y)
2

∞

dy + α

ze−α|x−y| un (y)dy

(5.8)

L

z2
+ z e−α|x−y| un+1 (y)dy + v n+1 (L)
˜
2

or

un+1 = α

∞

e−α|x−y|

−∞
∞

+α
L

e−α|x−y|

−

5 n
1
u (y) − un−1 (y) + un−2 (y) dy
4
4
z 2 3z
−
4
4

z
un+1 (y) + zun (y) − un−1 (y) dy + v n+1 (L).
˜
4

97

(5.9)

Here

v n+1 =
˜

−1
2α

δ(x − L)e−α|x−y| v n+1 (y)dy =

−1 −α|x−L| n+1
e
v
(L).
2α

Hence, due to the delta function, the term v is localized to the space-PML interface.
Consider the discretization of equation (5.8) or (5.9) using the midpoint rule. Let un
j
denote this fully discretized solution at (xj−1/2 , tn ), with xj−1/2 = −L + (j − 1/2)∆x.
The gridpoints corresponding to j = 1, 2, . . . N are in the computational region and those
corresponding to j = N + 1, N + 2, . . . N + M are in the PML region. The solution can then
be written as
N +M

un+1 =
j

Aji
i=1
N +m

un−1
un − i
i
2

Aji zun −
i

+
i=N +1

−α|xj−1/2 −xN | n+1
v

−e

z2
+ z un+1
i
2

(5.10)

or
N +M

un+1 =
j

Aji
i=1
N +m

+

1
5 n
−α|xj−1/2 −xN | n+1
ui − un−1 + un−2 − e
v
i
i
4
4

Aji
i=N +1

−

z 2 3z
−
4
4

z
un+1 + zun − un−1
i
i
4 i

(5.11)

with

n+1
−z un
n
v n+1 = e−z v n + λ un+1
N +1+p − uN −p
N +1+p − uN −p + e

98

.

The coefficient λ can be simplified using ν = α∆x, so that

σ ∆t
1
z
λ= 0
=
,
2α 2 (2p + 1)∆x
4(2p + 1)ν

and p is any non-negative integer of my choice to represent the derivative at x = L. The
matrix A is as given in section 2.3,

ν −|i−j|ν
Aji = 2 1 − e−ν/2 δji + 2 sinh
e
(1 − δji ),
2

1 ≤ j, i ≤ N + M

in terms of the Kronecker delta δji .
The important observation to be made here is that while the solution is given implicitly,
these implicit terms arise solely in the PML region. This will allow the continued use of the
fast summation methods used in the non-outflow algorithms, after further analysis of the
implicit equations which consitute the subsystem, N + 1 ≤ j ≤ N + M .

5.1.1.2

Implicit Update

The matrix A decomposes as

A = (b0 − b1 ) I + b1 T

where

b0 = 2 1 − e−ν/2 ,

99

b1 = 2 sinh

ν
2

and T is the Toeplitz matrix given in terms of d = e−ν by




. . . dN +M −1 


d
1
. . . dN +M −2 

,

.
.
.
...
.
.
.

.
.
.


dN +M −1 dN +M −2 . . .
1
1






T =





d

whose inverse is given by the tridiagonal matrix with perturbed corners

−d
1


2
1 −d 1 + d

−1 =
T

.
1 − d2  .
.
 .
.
.


.
.
0
.

...
...
..

.

..

.


0


0

.
.
.
.

1

Now, define

u = (uN +1 , . . . , uN +M )T ,
¯

then the subsystem corresponding to (5.10) or (5.11) can be rearranged and expressed by

I+

z2
1
+ z A un+1 = A (z + 1)¯n − un−1 − v n+1
¯
u
¯
¯
2
2

(5.12)

or

I+

z 2 3z
+
4
4

A un+1 = A
¯

z+

5
4

100

un −
¯

z
1
+ 1 un−1 + un−2
¯
¯
4
4

− v n+1
¯

with

vi = e
¯

−α xi−1/2 −xN

v n+1 ,

N + 1 ≤ i ≤ N + M.

We can rewrite the first order system as

I+

z2
1
+ z A un+1 = A (z + 1)¯n − un−1
¯
u
¯
2
2

− v n+1
¯

1
= ((1 + c0 )I + c1 T ) un+1 = (c0 I + c1 T ) (z + 1)¯n − un−1
¯
u
¯
2

− v n+1
¯

with

c0 = (b0 − b1 )

1 2
z +z ,
2

c1 = b1

1 2
z +z .
2

Multiply by T −1 to get

(1 + c0 )T −1 + c1 I un+1 = c0 T −1 + c1 I
¯

1
¯
(z + 1)¯n − un−1 − T −1 v n+1 .
u
¯
2

From here, note that the matrix on the left hand side, (1 + c0 )T −1 + c1 I is tri-diagonal
and thus this equation can be solved for un+1 using a tri-diagonal solver.
¯
An alternative approach is to bring all of the implicit terms to the left hand side of the
equation. We begin with the breakdown of the temporal derivative as

−α xi−1/2 −xN

v n+1 = e
¯

n+1
−z un
n
e−z v n + λ un+1
N +1+p − uN −p
N +1+p − uN −p + e

101

.

Rewrite this as

v n+1 = λED¯n+1 + e
¯
u

−α xi−1/2 −xN

e−z v n + λ e−z un +1+p − un −p
N
N

−α|xi−1/2 −xN |

where E is the matrix with e

,

on the ith diagonal entry and D is the matrix

with −1 in the (N − p)th column and 1 in the (N + p + 1)th column. Note that the matrix
multiplication ED results in a matrix with only two non-zero columns.
We can then rewrite (5.12) as

I+

z2
+ z A + λED un+1 =
¯
2

1
A (z + 1)¯n − un−1
u
¯
2

−e

−α xi−1/2 −xN

e−z v n + λ e−z un +1+p − un −p
N
N

.

Once again, reformulate A in terms of I and T as

1
¯
((1 + c0 )I + c1 T + λED) un+1 = (c0 I + c1 T ) (z + 1)¯n − un−1
¯
u
2
−e

−α xi−1/2 −xN

e−z v n + λ e−z un +1+p − un −p
N
N

Multiply by T −1 as before to get

1
(1 + c0 )T −1 + c1 I + λT −1 ED un+1 = (c0 T −1 + c1 I) (z + 1)¯n − un−1
¯
u
¯
2
− T −1 e

−α xi−1/2 −xN

e−z v n + λ e−z un +1+p − un −p
N
N

It must be noted here that, due to the λT −1 ED term, the left hand side is no longer
tri-diagonal. However, since T −1 is tri-diagonal and ED is a matrix with only two non-zero

102

columns, λT −1 ED will also be a tridiagonal matrix with two additional non-zero columns.
Thus, the left hand side will still be sparse. Furthermore, this has the benefit of creating
a right hand side that is truly explicit as it depends solely upon values from previous time
steps.
The second order system may be solved in a very similar way, resulting in either the
tri-diagonal system

I+

z 2 3z
+
4
4

A un+1 = A
¯

z+

5
4

un −
¯

z
1
+ 1 un−1 + un−2 − v n+1
¯
¯
¯
4
4

=⇒ ((1 + c0 )I + c1 T ) un+1 =
¯
(c0 I + c1 T )
=⇒

z+

5
4

un −
¯

z
1
+ 1 un−1 + un−2 − v n+1
¯
¯
¯
4
4

(1 + c0 )T −1 + c1 I un+1 =
¯
c0 T −1 + c1 I

z+

5
4

un −
¯

1
z
+ 1 un−1 + un−2
¯
¯
4
4

− T −1 v n+1 .
¯

Again, we have the option to remove all the implicit terms from the right hand side with

(1 + c0 )T −1 + c1 I + λT −1 ED un+1 =
¯
(c0 T −1 + c1 I)
− T −1 e

z+

5
4

−α xi−1/2 −xN

un −
¯

z
1
+ 1 un−1 + un−2
¯
¯
4
4

e−z v n + λ e−z un +1+p − un −p
N
N

.

Once again, the left hand side is now a tri-diagonal matrix with two additional non-zero
columns for every change in σ0 .

103

For the second order non-dissipative MOLT formulation, (5.6) becomes

1
(c∆t)2

un+1 − 2un + un−1 + σ∆t un+1 − un−1 + (σ∆t)2

un+1 + un−1
2

n−1
un+1 + uxx
= xx
− v n+1 .
2

We treat v in exactly the same fashion as above in the dissipative solvers. After reformulating this as a Helmholtz equation and inverting the operator, this leads to
∞

un+1 + un−1 = α

e−α|x−y| un (y)dy + α

−∞
∞

−α

∞

2ze−α|x−y| un−1 (y)dy

L

z 2 + z e−α|x−y| un+1 (y) + un−1 (y) dy + v n+1 (L),
˜

L

with

z=

σ∆t
.
2

This gives
N +M
n+1
uj

n−1
+ uj

Aji un − e
i

=

i=1
N +m

+

−α xj−1/2 −xN

Aji 2zun − z 2 + z
i

i=N +1

v n+1
un+1 + un−1
i
i

,

n+1
−2z un
n
v n+1 = e−2z v n + λ un+1
N +1+p − uN −p
N +1+p − uN −p + e

104

where

σ
∆t
z
λ= 0
=
.
α 2(2p + 1)∆x
(2p + 1)ν

The subsection corresponding to the PML region can then be expressed as

I + z2 + z A
vi = e
¯

un+1 + un−1 = A un + 2z un−1 − v ,
¯
¯
¯
¯
¯

−α xi−1/2 −xN

v n+1 ,

N + 1 ≤ i ≤ N + M.

Letting c0 = (b0 − b1 ) z 2 + z and c1 = b1 z 2 + z , this now becomes

((1 + c0 ) I + c1 T ) un+1 + un−1 = (c0 I + c1 T ) un + 2z un−1 − v .
¯
¯
¯
¯
¯

Multiply both sides by T −1 to get the tri-diagonal system

(1 + c0 ) T −1 + c1 I

un+1 + un−1 = c0 T −1 + c1 I
¯
¯

un + 2z un−1 − T −1 v .
¯
¯
¯

Furthermore, we can obtain the system with only explicit terms on the right hand side
as

I + z 2 + z A + λED
A un + 2z un−1 + e
¯
¯

un+1 + un−1 − λED¯n−1 =
¯
¯
u

−α xi−1/2 −xN

e−2z v n + λ e−2z un +1+p − un −p
N
N

105

.

This implies that

((1 + c0 ) I + c1 T + λED) un+1 + un−1 = λED¯n−1 + (c0 I + c1 T ) un + 2z un−1
¯
¯
u
¯
¯
−α xi−1/2 −xN

+e

e−2z v n + λ e−2z un +1+p − un −p
N
N

.

Multiplication by T −1 then gives

(1 + c0 ) T −1 + c1 I + λT −1 ED

un+1 + un−1 =
¯
¯

λT −1 ED¯n−1 + c0 T −1 + c1 I
u
+ T −1 e

5.1.1.3

−α xi−1/2 −xN

un + 2z un−1
¯
¯

e−2z v n + λ e−2z un +1+p − un −p
N
N

.

Performance

While using a constant damping factor σ(x) = σ0 H(x − L) may be fairly easy to implement,
it does not necessarily provide a good method for truncating the domain. The biggest issue
is that, due to the Heaviside function, σ(x) is not smooth when switching from the interior
of the domain to the PML region and thus, the regions are no longer ‘perfectly’ matched.
This break in continuity may cause reflections of the wave as it passes the boundary into
the absorbing layer. To check for these reflections and to quantify the damping potential of
a piece-wise constant damping function, we start with an initial Gaussian pulse and let it
propagate for a full pass through the domain with a perfectly matched region on the right
side and see what percentage of the initial function still remains. The parameters were:
2

• Initial conditions are f (x) = e−100(x−0.5) , g(x) = 0.
• The mesh sizes are ∆x = ∆t = 0.1 with a wave speed of c = 1.
106

• The interior domain is given by [−1, 1].
• The damping layer is [1, 5].
• The initial damping constant was given by one of {0.1, 0.01, 0.001, 0.0001}.
• The damping constant was doubled halfway through damping region.
At time t = 1, a check is perfomed on both ||u(x)||2 and ||u(x)||∞ in the range [0, 1].
At this point, it is safe to assume that the right bound wave has entered the PML region
while the left bound wave has not yet reached the matched boundary. This will allow us to
compare the spurious numerical reflection coming off of the boundary from the right bound
wave caused by the discontinuity in the damping function. Finally, at a time of t = 12, which
corresponds to the time necessary for the wave to make one full pass through the domain, we
once again check both ||u(x)||2 and ||u(x)||∞ to see how much of the initial pulse remains.
The results are shown in tables 5.1 and 5.2.
As can easily be seen from the tables, the reflection caused by the discontinuity in the
damping function may eliminate much of the benefit that comes from the actual damping.
In addition, the actual damping itself does not perform well, removing no more than 70% of
the wave despite the damping region taking up well over half of the domain.

107

Table 5.1: Reflection and Damping for the Dissipative 1st Order Scheme with a Piecewise
Constant Damping Function
σ0
0
0.0001
0.001
0.01
0.1

∞-norm at t = 1
0.000014176
0.000016310
0.000035517
0.00022812
0.0022091

2-norm at t = 1
0.000022378
0.000029892
0.00018632
0.0018358
0.018325

∞-norm at t = 12
0.17580
0.17559
0.17374
0.15625
0.054052

2-norm at t = 12
1.48488
1.48313
1.46747
1.31976
0.45717

Table 5.2: Reflection and Damping for the Dissipative 2nd Order Scheme with a Piecewise
Constant Damping Function
σ0
0
0.0001
0.001
0.01
0.1

5.1.2

∞-norm at t = 1
0.00000098983
0.0000039184
0.000030283
0.00029472
0.0030214

2-norm at t = 1
0.0000013516
0.000020798
0.00020579
0.0020593
0.020930

∞-norm at t = 12
0.35042
0.35001
0.34638
0.31205
0.11020

2-norm at t = 12
2.44360
2.44077
2.41540
2.17588
0.76839

Polynomial Damping Functions

To correct the matching deficiencies of the piece-wise constant damping functions from the
previous section, consider damping functions of the form

σ(x) = σ0

x−L m
,
d

m ≥ 2.

We make use of a Heaviside function to extend this to the entire domain as

σ(x) = σ0

x−L m
H(x − L),
d

m ≥ 2.

These polynomial damping functions will ensure that σ(x) ∈ Cm .

108

(5.13)

5.1.2.1

Implementation

We now wish to create a fast implementation of equations (5.4) and (5.5) incorporating
(5.13). We first note that by use of an integrating factor, the solution to (5.5) is given by
t

v(x, t) = σ (x)
=⇒

v(x, tn+1 )

=⇒

v(x, tn+1 )

= σ (x)

e−σ(x)(t−s) ux (x, s)ds

0

tn+1 −σ(x) tn+1 −s

e

0

(m+1)∆t

n

=

σ (x)

ux (x, s)ds

e−σ(x)((n+1)∆t−s) ux (x, s)ds.

m∆t

m=0

We then approximate each integral using trapezoidal quadrature and arrive at

v n+1

σ (x)∆t
=
2

n

ux (x, m∆t)eσ(x)(m−n−1)∆t

(5.14)

m=0

+ ux (x, (m + 1)∆t)eσ(x)(m−n)∆t .

Applying the standard first order discretization to equation (5.4) gives
1 + 2σ(x)∆t + σ(x)2 ∆t2 un+1 (x)
∆t2
=

c2 un+1 (x) − σ
xx

tn+1

(x)

0

−

2(1 + ∆tσ(x))un
∆t2

+

un−1 (x)
∆t2

(5.15)

ux (x)e−σ(x)(t−s) .

Combining equations (5.15) and (5.14) produces the following scheme which we would

109

like to solve quickly:
1 + 2σ(x)∆t + σ(x)2 ∆t2 un+1 (x)
∆t2
σ (x)∆t
−
2

−

2(1 + ∆tσ(x))un
∆t2

+

un−1 (x)
∆t2

= c2 un+1 (x)
xx

n

ux (x, m∆t)eσ(x)(m−n−1)∆t + ux (x, (m + 1)∆t)eσ(x)(m−n)∆t .
m=0

This then may be re-arranged into a Helmholtz equation as

∂xx − α2 un+1 = α2

2σ∆t + σ 2 ∆t2 un+1 − 2(1 + σ∆t)un + un−1 + v n+1 .

An inversion of the Helmholtz operation then yields

un+1 = −

α
2

e−α|x−y|

2σ(y)∆t + σ 2 (y)∆t2 un+1 (y)

(5.16)

− 2(1 + σ(y)∆t)un (y) + un−1 (y) dy + v n+1
˜

with

v n+1 = −
˜

1
2α

e−α|x−y|

σ (y)∆t
2

(5.17)

n

uy (y, m∆t)eσ(y)(m−n−1)∆t + uy (y, (m + 1)∆t)eσ(y)(m−n)∆t dy.

·
m=0

Define

ω0 (y) =

σ (y)∆t −α|x−y|
e
,
4α

110

ωk+1 = e−σ(y)∆t ωk (y)

such that integration by parts on (5.17) will yield

v n+1
˜

1
=−
4α
n

n

∞

m=0 L
∞

σ (y)∆te−α|x−y|−(n−m)∆tσ(y) um+1 (y) + e−∆tσ(y) um (y) dy
y
y

ωn−m (y)um+1 (y) + ωn+1−m (y)um (y)dy
y
y

=−
m=0 L
n

ωn−m (y)um+1 (y) + ωn+1−m (y)um (y)

=−
m=0
n

∞

+
m=0 L

∞
L

ωn−m (y)um+1 (y) + ωn+1−m (y)um (y)dy.

The left boundary term will vanish due to the stipulation that the minimum order of our
polynomial was at least 2, which guarantees ωk (L) = 0, k ∈ {0, 1, . . . , n}. Additionally, the
right boundary term will vanish due to the negative exponential. This leaves that
n

v n+1
˜

∞

=
m=0 L

ωn−m (y)um+1 (y) + ωn+1−m (y)um (y)dy.

Furthermore, for a truncated domain, we can ensure the right boundary term disappears
by enforcing zero Dirichlet boundary conditions at the end of the damping region.
Using the same techniques as in section 5.1.1.1, we now consider the fully discretized
form of equation (5.16), given by
N +M

un+1 (x

j)

=

1
Aji un − un−1
i
2 i

i=1
N +m

−
i=N +1

(2) 1
(1)
(1)
Σji un+1 + Σji un+1 − Σji un
i
i
2

111

(5.18)
+ vj ,
˜n+1

where

(1)

Σji = α
(2)

Σji = α

xi

σ(y)∆te

xi−1
xi

−α xj−1/2 −y

(σ(y)∆t)2 e

dy and

−α xj−1/2 −y

(5.19)

dy,

xi−1

which are defined for 1 ≤ i, j ≤ N + M , but noting that σ(y) = 0 outside of the PML region.
Further, we can simplify the term vj . After application of the midpoint rule on the u
˜n+1
terms, we have
N +M

n

vj
˜n+1 =

L+i∆x

ωn−m (y)um+1 + ωn+1−m (y)um dy
i
i

m=0 i=N +1 L+(i−1)∆x
n N +M
x
[ωn−m (y)]xi um+1
i−1 i
m=0 i=N +1

=

x

+ [ωn+1−m (y)]xi um .
i−1 i

This implies that
N +M

vj
˜n+1 =

n

x

σ (y)∆t −α xj−1/2 −y −(n−m)σ(y)∆t i
e
e
um+1
i
4α
xi−1
i=N +1 m=0
N +M

n

x

σ (y)∆t −α xj−1/2 −y −(n+1−m)σ(y)∆t i
e
e
um
+
i
4α
xi−1
i=N +1 m=0
N +M

=

·

σ (xi )∆t −α xj−1/2 −xi
e
4α

i=N +1
n
e−(n−m)σ(xi )∆t

um+1 (xi ) + e−σ(xi )∆t um (xi )

m=0

−

σ (xi−1 )∆t −α xj−1/2 −xi−1
e
4α

n

·

e−(n−m)σ(xi−1 )∆t um+1 (xi ) + e−σ(xi−1 )∆t um (xi ) .

m=0

112

Therefore,
N +M

vj
˜n+1

=

−α xj−1/2 −xi

e
i=N +1
N +M

=
i=N +1

n+1
v1,i − e

−α xj−1/2 −xi−1

n+1
v2,i

(1) n+1
(2) n+1
Bji v1,i − Bji v2,i

where

(1)

Bji = e−ν|j−i−1/2| ,

(2)

Bji = e−ν|j−i+1/2| ,

and the sums are defined by the auxiliary variables

n+1
v1,i

σ (xi )∆t
=
4α

n+1
v2,i =

n

um+1 (xi ) + e−σ(xi )∆t um (xi ) ,

m=0
n

σ (xi−1 )∆t
4α

um+1 (xi ) + e−σ(xi−1 )∆t um (xi ) .

m=0

We can simplify this further by using a recurrence relation to remove the sum. Thus, the
auxiliary variables can be updated locally in time by

n+1
n
v1,i = Ei v1,i + Fi un+1 + Ei un ,
i
i

n+1
n
v2,i = Ei−1 v2,i + Fi−1 un+1 + Ei−1 un
i
i

(5.20)

where

Ei = e−σ(xi )∆t ,

Fi =

σ (xi )∆t
,
4α

N ≤ i ≤ N + M.

(5.21)

Thus, for each time step, we update the equations (5.20) as well as solving for equation

113

(5.18) in the following form
N +M

un+1 (x

j)

=

Aji
i=1
N +m

un
i

1 n−1
− ui
2

N +m

−
i=N +1

(1)
(2) 1
(1)
Σji un+1 + Σji un+1 − Σji un
i
i
2

(1)

+
i=N +1
N +m

−
i=N +1

n
Ei v1,i + Fi un+1 + Ei un
i
i

(2)

n
Ei−1 v2,i + Fi−1 un+1 + Ei−1 un
i
i

(5.22)

Bji
Bji

.

We now specifically consider the damping function

σ(x) = σ0

x−L 2
d

for the purpose of fast implementation of this scheme. The techniques that will follow can
be adapted in a very straight-forward manner for polynomials of higher order. However, we
have chosen to use a polynomial of degree 2 as it will be the easiest to implement as well as
performing the best as a damping function as can be seen in figure 5.1.
Since we have assumed the damping region to be M grid points wide, we have d = M ∆x
and can write

si = σ(xi )∆t = σ0 ∆t

i−N 2
,
M

N ≤ i ≤ N + M.

Similarly,

si = σ (xi )∆t = 2σ0 ∆t

i−N
M

114

,

N ≤ i ≤ N + M.

Damping Profiles for Low Order Polynomials
1

0.8

Order 2
Order 3
Order 4

0.6

0.4

0.2

0
0

0.1

0.2

0.3

0.4

0.5
(x-L)/d

0.6

0.7

0.8

0.9

Figure 5.1: Damping Profiles for Low Order Polynomials
From these expressions, we can now evaluate (5.19)
xi

(1)

Σji = ασ0 ∆t

xi−1
1/2

= νσ0 ∆t

−1/2

y − L 2 −α xj−1/2 −y
e
dy
d
i − 1/2 − N + z 2 −ν|j−i−z|
e
dz
M

using the change of variables z = y − xi−1/2 /dx. This then becomes
σ ∆t 1/2
(1)
Σji = ν 0
z 2 + (i − N − 1/2)2 + (2i − 2N − 1)z e−ν|j−i−z| dz.
2 −1/2
M
Note that Σ(1) can be decomposed into terms of the form

(k)
Gji

1/2

=ν

z k e−ν|j−i−z| dz.

−1/2

115

1

We now evaluate these expressions for each value of k and the cases i = j and i = j.
When k = 0, we recover

(0)

Gji = Aji = 2 1 − e−ν/2 δji + 2 sinh

ν −|i−j|ν
e
(1 − δji ).
2

More generally, for the diagonal terms and arbitrary k,
1/2

(k)

Gii = ν

=

=

z k e−ν|z| dz

−1/2




0,


k

 2 ν/2



wk e−w dw, k
k 0
ν


0,

 2k!

ν −ν/2


e
1 − bk
k
2
ν

odd
even
k

odd

k

even,

where
k

bk (z) =
=0

ν
Note that bk
2

z
!

.

(k)
is exactly the Taylor polynomial for eν/2 of order k, so that Gii is

well defined as ν tends to 0.

116

For j = i, we similarly have

(k)

Gji = νe−|j−i|ν
= e−|j−i|ν

1/2

z k e−νz dz

1

−1/2
ν/2

νk

−ν/2

wk e− sgn (j−i)w dw

ν −ν/2
ν ν/2
e
− bk
e
2
2
νk
2k!
ν
ν
ν
ν
= e−|j−i|ν (sgn (j − i))k
ck
sinh
− sk
cosh
k
2
2
2
2
ν
= e−|j−i|ν

k!

sgn(j − i)k bk −

,

where ck (z) and sk (z) are the even and odd terms, respectively, of bk (z). Because of the
sign change due to odd values of k, it will be useful to define the upper and lower triangular
Toeplitz matrices U and L = U T , where


d2

...
 0 d


...

0 d



... ...
U =




0



dN +M −1
.
.
.
.
.
.
d
0

















with d = e−ν . Further, we have the Toeplitz matrix







−1 = I + U + L = 
T










. . . . . . dN +M −1 


.
.. ..
.

.
.
d
1
.



.
.
.. .. ..
.
.

.
.
.
.
.


.

.. ..
.
.
. 1

.
d


dN +M −1 . . . . . . d
1
1

117

d

which is the inverse of the tridiagonal matrix




−d
1





..
−d 1 + d2

.




1 

..
..
..
T =

.
.
.
.

1 − d2 




..
. 1 + d2 −d





−d
1
Assembling all of this, we can now write

even



g (U − L),
 k1

G(k) =




g I + g (U + L), k
 k0
k1

odd

k

where

gk0 =

2k!
νk

1 − bk

ν −ν/2
e
,
2

gk1 =

2k!
νk

ck

ν
ν
ν
ν
sinh
− sk
cosh
2
2
2
2

.

Having fully determined expressions for G(k) , we now express

σ ∆t
Σ(1) = 0
G(2) + 2G(1) D + G(0) D2
M2
σ ∆t
= 0
g20 I + g00 D2 + U g21 I + 2g11 D + g01 D2 + L g21 I − 2g11 D + g01 D2
2
M
= T −1 Σ1 + U Σ1 − Σ1 + L Σ1 − Σ1
0
1
0
2
0

(5.23)

where we have used the fact that I = T −1 − U − L and D is the diagonal matrix with

118

elements

1
1 3
. Additionally, we have defined
0, . . . , 0, , , . . . , M −
2 2
2
σ ∆t
Σ1 = 0
g20 I + g00 D2 ,
0
2
M
σ ∆t
g21 I + 2g11 D + g01 D2 ,
Σ1 = 0
1
2
M
σ ∆t
Σ1 = 0
g21 I − 2g11 D + g01 D2 .
2
2
M

Likewise,

Σ(2) =

σ0 ∆t 2

=

σ0 ∆t 2

M2
M2

G(4) + 4G(3) D + 6G(2) D2 + 4G(1) D3 + G(0) D4
g40 I + 6g20 D2 + g00 D4

+ U g41 I + 4g31 D + 6g21 D2 + 4g11 D3 + g01 D4
+ L g41 I − 4g31 D + 6g21 D2 − 4g11 D3 + g01 D4
= T −1 Σ2 + U Σ2 − Σ2 + L Σ2 − Σ2
0
1
0
2
0

,

with

Σ2 =
0
Σ2
1

=

Σ2 =
2

σ0 ∆t 2
M2

σ0 ∆t 2
M2
σ0 ∆t 2
M2

g40 I + 6g20 D2 + g00 D4 ,
g41 I + 4g31 D + 6g21 D2 + 4g11 D3 + g01 D4 ,
g41 I − 4g31 D + 6g21 D2 − 4g11 D3 + g01 D4 .

119

(5.24)

We can also write the matrices B (1) and B (2) using combinations of T −1 , U and L as

ν
L,
2
ν
B (2) = e−ν/2 (I + L) + eν/2 U = e−ν/2 T −1 + 2 sinh
U.
2

B (1) = e−ν2 (I + U ) + eν/2 L = e−ν/2 T −1 + 2 sinh

(5.25)

Due to the symmetry and structure of T ,




0 d





2 ...
0 −d





1 

..
..
TU =

,
.
.
2

1−d 



2

−d
d 




0 −d2
and similarly




−d2




 d

1 

TL =

2
1−d 





0




2 ...

−d



... ...
.


2 0

−d


d 0

This makes it possible to write expressions for T B (1) , etc. as a combination of diagonal
matrices and products with T U and T L. From this, we will show that equation (5.22) can
be written as a tridiagonal system.
Upon collecting common terms in equation (5.22), and considering only N + 1 ≤ j ≤

120

N + M , we have

1
I + Σ(1) + Σ(2) − B (1) F + + B (2) F − un+1
2

(5.26)

1
= A + Σ(1) + B (1) F + E + − B (2) F − E − un − Aun−1
2
n
n
+ B (1) E + v1 − B (2) E − v2 ,

where the plus and minus superscripts of E and F indicate whether all but the first or last
entry of the M + 1 diagonal entries given by equations (5.21) is taken. Multiply by T and
make use of the decompositions given by (5.23), (5.24) and (5.25), so that

n
n
T2 un+1 = T1 un − T0 un−1 + T+ v1 − T− v2 ,

where

T2 = T

1
I + Σ(1) + Σ(2) − B (1) F + + B (2) F −
2

1
= T + Σ1 + Σ2 − e−ν/2 (F + − F − )
0
2 0
1
1
ν
− T U Σ1 − Σ1 + Σ2 − Σ2 − 2 sinh
F−
0
1
0
1
2
2
2
1
1
ν
− T L Σ1 − Σ1 + Σ2 − Σ2 − 2 sinh
F+ ,
0
2
0
2
2
2
2

121

T1 = T A + Σ(1) + B (1) F + E + − B (2) F − E −
= T A + Σ1 + e−ν/2 (E + F + − E − F − )
0
ν
− T U Σ1 − Σ1 − 2 sinh
F +E +
1
0
2
ν
F −E − ,
− T L Σ1 − Σ1 − 2 sinh
2
0
2

T0 =

TA
,
2

T+ = e−ν/2 E + + 2 sinh

ν
T L E+ ,
2

T− = e−ν/2 E − + 2 sinh

ν
T U E− .
2

and

The system is solved with a basic tridiagonal algorithm for un+1 in the PML region. It
j
n+1 and v n+1 .
is then used to update v1
2

In the free space region, note that equation (5.22) gives the solution for un+1 explicitly
j
for 1 ≤ j ≤ N , since the sums involving PML terms will not be carried over this range.
Thus, all terms arising from the PML will only appear as a constant vector in the N × N
subsystem. Furthermore, this vector has a simple form due again to the structure of T , U
and L.
If we instead choose to use a second order dissipative discretization for equation (5.4),

122

we now recover
2 + 3σ(x)∆t + σ(x)2 ∆t2 un+1 (x)

−

∆t2
+

(4 + σ(x)∆t)un−1 (x)
∆t2

−

un−2 (x)
∆t2

=

(5 + 4σ(x)∆t)un
∆t2
c2 un+1 (x) − σ
xx

tn+1

(x)

0

ux (x)e−σ(x)(t−s)

instead, but still treat

v n+1

tn+1

= σ (x)

0

−σ(x) tn+1 −s

ux (x)e

in exactly the same manner. Once again, rearrange this into a Helmholtz equation as

∂xx − α2 un+1 =

α2
2

3σ∆t + σ 2 ∆t2 un+1 − (5 + 4σ∆t) un

+ (4 + σ(x)∆t) un−1 − un−2 + v n+1 ,
√

with α =

2
here. Invert the Helmholtz operator to get
c∆t

un+1 (x) = −

α
4

e−α|x−y|

3σ(y)∆t + σ 2 (y)∆t2 un+1 (y)

− (5 + 4σ(y)∆t)un (y) + (4 + σ(y)∆t)un−1 (y) − un−2 (y) dy + v n+1 .
˜

with v n+1 defined as before.
˜

123

Discretization using the midpoint approximation will then result in
N +M

un+1 (x

j)

=

Aji
i=1

N +m

−
i=N +1

5 n
1
ui − un−1 + un−2
i
4
4 i

3 (1) n+1 1 (2) n+1
1 (1)
(1)
Σji ui
˜n+1
+ Σji u
− Σji un + Σji un−1 + vj ,
i
4
4
4

with Σ(1) and Σ(2) defined as before.
Following the exact same procedure as for the 1st order dissipative scheme will yield that
equation (5.26) now becomes

1
3
I + Σ(1) + Σ(2) − B (1) F + + B (2) F − un+1
4
4
=

5
A + Σ(1) + B (1) F + E + − B (2) F − E − un
4

1
1
n
n
− A + Σ(1) un−1 + Aun−2 + B (1) E + v1 − B (2) E − v2 .
4
4

Multiply by T to give that all the matrices on the left hand side become tri-diagonal at
worst, which gives

n
n
T3 un+1 = T2 un − T1 un−1 + T0 un−2 + T+ v1 − T− v2 ,

124

where

3
1
I + Σ(1) + Σ(2) − B (1) F + + B (2) F −
4
4

T3 = T

1
3
= T + Σ1 + Σ2 − e−ν/2 (F + − F − )
0
4
4 0
− TU

3 1 3 1 1 2 1 2
ν
Σ0 − Σ1 + Σ0 − Σ1 − 2 sinh
F−
4
4
4
4
2

− TL

ν
3 1 3 1 1 2 1 2
Σ0 − Σ2 + Σ0 − Σ2 − 2 sinh
F+ ,
4
4
4
4
2

5
A + Σ(1) + B (1) F + E + − B (2) F − E −
4

T2 = T
=

5
T A + Σ1 + e−ν/2 (E + F + − E − F − )
0
4

ν
F +E +
2
ν
− T L Σ1 − Σ1 − 2 sinh
F −E − ,
0
2
2
− T U Σ1 − Σ1 − 2 sinh
0
1

T1 = T

1
A + Σ(1)
4

= TA +

1
Σ1 − T U Σ1 − Σ1 − T L Σ1 − Σ1
0
0
1
0
1
4

1
T0 = T A,
4

T+ = e−ν/2 E + + 2 sinh

125

ν
T L E+ ,
2

,

and

T− = e−ν/2 E − + 2 sinh

ν
T U E− .
2

We can then find un+1 in the PML region using a tri-diagonal solver before explicitly
updating for un+1 in the entire domain.
For the second order non-dissipative scheme, after the standard discretization, equations
(5.4) and (5.5) become
1 + 2s + 2s2 w − 2un − 4sun−1
(c∆t)2

1
= wxx − v n+1 ,
2

n

v n+1

e−2s(n−m) um+1 + e−2s um ,
x
x

=s
m=0

where

w = un+1 + un−1 ,

s(x) =

σ(x)∆t
.
2

Arranging this in a Helmholtz type equation gives

∂xx − α2 w = −2α2 un + 2sun−1 − s + s2 w + v n+1
√

with α =

2
. After inverting the operator, this now yields
c∆t

w(x) = α

e−α|x−y| un (y) + 2s(y)un−1 (y) − s(y) + s2 (y) w(y) dy + v n+1
˜

126

with

v n+1 (x)
˜

1
=−
2α

n

e−α|x−y| s

e−2(n−m)s(y) um+1 (y) + e−2s(y) um (y) dy.
y
y
m=0

After a full discretization of the domain, we obtain
N +m

N +M

Aji un
i

wj =

+
i=N +1

i=1

(1)
(1)
(2)
2Σji un−1 − Σji + Σji wi + vj .
˜n+1
i

Handling the temporal integral vj
˜n+1 in the same way as for the dissipative cases then
gives
N +M

N +m

Aji un
i

wj =
i=1
N +m

+
i=N +1
N +m

−
i=N +1

+
i=N +1

(1)
(2)
(1)
2Σji un − Σji + Σji wi
i

(1)

n
Ei v1,i + Fi Ei un + wi − un−1
i
i

(2)

n
Ei−1 v2,i + Fi−1 Ei−1 un + wi − un−1
i
i

Bji

Bji

.

This then becomes

I + Σ(1) + Σ(2) − B (1) F + + B (2) F − w
= A + B (1) E + F + − B (2) E − F − un + 2Σ(1) − B (1) F + + B (2) F − un−1
n
n
+ B (1) E + v1 − B (2) E − v2 .

127

Multiplication by T then produces

n
n
Tw w = T1 un + T0 un−1 + T+ v1 − T− v2 + S,

where

Tw = I + Σ1 + Σ2 − e−ν/2 (F + − F − )
0
0
− T U I + Σ1 − Σ1 + Σ2 − Σ2 − a1 F −
1
0
1
0
− T L I + Σ1 − Σ1 + Σ2 − Σ2 + a1 F + ,
0
2
0
2

T1 = a0 I + e−ν/2 (E + F + − E − F − )
− T U (a0 − a1 )I + a1 E − F −
− T L (a0 − a1 )I − a1 E + F + ,

T0 = 2Σ1 + 2Σ2 − e−ν/2 (F + − F − )
0
0
− T U 2Σ1 − 2Σ1 + 2Σ2 − 2Σ2 − a1 F −
0
1
0
1
− T L 2Σ1 − 2Σ1 + 2Σ2 − 2Σ2 + a1 F + ,
0
2
0
2

T+ = e−ν/2 E + + a1 T L E + ,

128

(5.27)

T− = e−ν/2 E − + a1 T U E − ,

and finally the vector S is nonzero only in the first entry, and is comprised of the matrix
vector product over j = 1, 2, . . . N
N

(T L)N +1,i un ,
i

Sj = δj,N +1

N + 1 ≤ j ≤ N + M.

i=1

Equation (5.27) may then be solved in the PML region using a tri-diagonal solver. This
solution may then be used to calculate an explicit solution for un+1 over the entire domain.

5.1.2.2

Optimization

Optimizing the damping of a wave in a PML layer has been a subject of significant study, but
most authors believe that the result is very much problem specific. Thus, we will conduct a
numerical study for optimizing the absorbance of a perfectly matched layer for a quadratic
damping function. The variables will consist of both σ0 and the length of the PML region.
The other parameters of the study are as such:
2

• The initial conditions are f (x) = e−100x and g(x) = 0.
• The interior domain is [−1, 1].
• The PML region has a length of {0.5, 2, 3}.
• ∆x = ∆t = 0.01 with a wave speed of c = 1.
• The chosen final time is tf = 20, such that all waves would end on a recombination.
• σ0 = {1, 2, 3, 4, 5, 10, 15, 20, 25}.
129

Table 5.3: Optimization of the Dissipative 1st Order Quadratic Damping Function

σ0
0
1
2
3
4
5
10
15
20
25

d = 0.5
∞-norm
2-norm
0.14901
1.35786
0.046120
0.40325
0.014676
0.12161
0.0046670
0.036095
0.0014033
0.010184
0.00035072 0.0031107
0.000045637 0.00032578
0.000044945 0.00032420
0.000038350 0.00027405
0.000030196 0.00023973

d=2
∞-norm
2-norm
0.13877
1.30731
0.021976
0.19768
0.0047455 0.045587
0.0014279 0.018349
0.00086525 0.010557
0.00045273 0.0058206
0.00044261 0.0056118
0.00047148 0.0044793
0.00053009 0.0050003
0.00038082 0.0032459

d=3
∞-norm
2-norm
0.13422
1.28520
0.016763
0.16238
0.0029356 0.035231
0.00087016 0.011461
0.0010083 0.011400
0.0010008 0.011423
0.00065305 0.0076492
0.00050443 0.0053598
0.00050236 0.0066065
0.00036747 0.0044241

Both the 2-norm and ∞-norm of u was taken at the end of the run. A baseline with σ0 = 0
was run for each PML length to determine how much damping had actually occurred. The
results can be seen in table 5.3. From the table, while it appears that the first few increases of
σ0 are quite helpful, the results very much level out shortly thereafter. This could be due to
the larger co-efficient exacerbating the error that is caused by the numerical approximation
of the integral. The best scenario appears to be taking a shorter PML region accompanied
with a relatively larger σ0 , especially in consideration of the fact that a shorter PML region
will constitute less work being done by the tri-diagonal solver when solving the implicit step
which appears only in the PML region.

5.2

One-Way Waves

An alternative approach to the idea of a perfectly matched damping layer is to construct
a boundary that appears invisible that will simply let the outward propagating wave pass
through without creating any reflection. While this may be rather difficult to do in the
general case, it is quite feasible in the one dimensional case under a few basic assumptions.
130

Assume that the initials conditions f (x) and g(x) in equation (1.1) are of compact support
and that this support does not extend beyond Ω0 , the desired domain of computation. For
the sake of simplicity, let Ω0 = [−L, L].
The first order dissipative MOLT solution to the 1D wave equation is then given by

1
α un (y) − un−1 (y) e−α|x−y| dy,
2
Ωt n

un+1 (x) =

(5.28)

where Ωtn is the domain of dependence for u(x) at time tn .
Since the wave speed propagation is limited by c, it is easy to see that Ωtn is fully
contained in [−L − ctn , L + ctn ], which gives that (5.28) now becomes

un+1 (x) =

L+ctn

1
α un (y) − un−1 (y) e−α|x−y| dy.
2
−L−ctn

(5.29)

This is equivalent to the free space solution as un+1 (x) will be identically zero outside of
this region.
We now wish to compute the solution only over the initial region Ω0 while still factoring in
the contribution from the additional domain of Ωtn+1 /Ω0 . This is accomplished by writing
(5.29) as

un+1 (x) =
n+1
S1

L

1
n+1
n+1
α un (y) − un−1 (y) e−α|x−y| dy + S1 + S2 ,
2
−L

(5.30)

−L

=

n+1 =
S2

1
α un (y) − un−1 (y) e−α(x−y) dy,
2
−L−ctn
L+ctn
L

1
α un (y) − un−1 (y) e−α(y−x) dy.
2

Since we are considering the case of outflow boundary conditions, the key step here is to
131

restrict only to outgoing waves that will not reflect back into the domain. This means that
we must impose non-reflecting boundary conditions at x = ±L. For this one dimensional
case, this condition simplifies to the respective transport equations

ut + cux = 0,

x≥L

ut − cux = 0,

x ≤ −L.

We now consider the case where x ≥ L and the case where x ≤ −L will be similar.
Tracing along the characteristics gives that

u(L + y, t) = u L, t −

y
c

for y > 0.

Thus, for x < L, we can write

n+1 =
S2

L+ctn
L
L+ctn

=
L

1
α un (y) − un−1 (y) e−α(y−x) dy
2
1
α u(y, tn ) − u(y, tn−1 ) e−α(y−x) dy
2

= αe−α(L−x)
= αce−α(L−x)

ctn
0
tn
0

1
u(L + y, tn ) − u(L + y, tn−1 ) e−αy dy
2
1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds.
2

This, in effect, has transformed the spatial integral into a temporal integral and therefore,
we must only track the time history at the boundary points x = ±L. Additionally, we may
proceed to track this temporal integral using a recurrence relation which limits storage to

132

the previous two time steps. Note that
tn

n+1 = e−α(L−x) αc
S2

1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds
2

0

˜n+1
= e−α(L−x) S2 .

Making use of the change of variables s = s + ∆t gives

˜n+1 = αc
S2

tn
0

1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds
2

˜n
= e−αc∆t S2 + αc
˜n
= e−1 S2 + αc

∆t
0

∆t
0

1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds
2

1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds.
2

∆t
We now choose to integrate 0 u (L, tn − s) e−αcs ds by approximating u(L, tn − s) with

a linear function through points at u(L, tn ) and u(L, tn−1 ). This linear approximation is
given by

u(L, tn − s) ≈

u(L, tn−1 ) − u(L, tn )
s + u(L, tn ) = zs + b.
∆t

133

(5.31)

Thus,
∆t

αc

0

s
∆t
−
1
u (L, tn
=
(zs + b)e ∆t ds
∆t 0
s
s ∆t
∆t
−
−
1
z∆te ∆t ds
=
−(zs + b)∆te ∆t
+
∆t
0
0
− s) e−αcs ds

1
=
∆t
=

s ∆t
s ∆t
2 ze− ∆t
−(zs + b)∆te ∆t
− ∆t
−

0

0

1
−(z∆t + b)∆te−1 − ∆t2 ze−1 + b∆t + z∆t2
∆t

= −(z∆t + b)e−1 − ∆tze−1 + b + z∆t
= − un−1 (L) − un (L) + un (L) e−1 − un−1 (L) − un (L) e−1
+ un (L) + un−1 (L) − un (L)
= e−1 un (L) + (1 − 2e−1 )un−1 (L).

This also gives that
∆t

αc

0

u (L, tn−1 − s) e−αcs ds = e−1 un−1 (L) + 1 − 2e−1 un−2 (L).

134

Thus, estimating this integral using the trapezoidal rule results in
∆t

αc

0

1
u (L, tn − s) − u (L, tn−1 − s) e−αcs ds
2

1 −1 n−1
e u
(L) + 1 − 2e−1 un−2 (L)
2
1 − 2e−1 n−2
5 −1
un−1 (L) −
= e−1 un (L) + 1 − e
u
(L)
2
2
= e−1 un (L) + 1 − 2e−1 un−1 (L) −

=

αc∆t
2

=

1
2

1
1
e−αc∆t u(L, tn−1 ) − u(L, tn−2 ) + u(L, tn ) − u(L, tn−1 )
2
2

u(L, tn ) + e−1 −

1
2

u(L, tn−1 ) −

e−1
u(L, tn−2 ) .
2

This implies that
−1
˜2 = e−1 S 2 + e−1 un (L) + 1 − 5 e−1 un−1 (L) − 1 − 2e un−2 (L).
˜
Sn+1
n
2
2

Similarly,

n+1
˜n+1
S1 (x) = e−α(L+x) S1

˜n+1 = e−1 S n + e−1 un (−L) + 1 − 5 e−1 un−1 (−L)
˜
S1
1
2
−

1
1 − 2e−1 n−2
u
(−L)
2
2

u(−L, tn ) + e−1 −

135

1
2

u(−L, tn−1 ) −

e−1
u(−L, tn−2 ) .
2

Thus, the fully discretized solution over Ω0 will be given by
m

un+1
j

=
i=1

m

=
i=1
m

=
i=1

1 n−1
Aji un − ui
i
2

n+1
n+1
+ S1 + S2

1 n−1
−α(L+xj ) ˜n+1
−α(L−xj ) ˜n+1
+e
S1 + e
S2
Aji un − ui
i
2
1 n−1
Aji un − ui
i
2

+e

−α(L+xj )

5
1 − 2e−1 n−2
˜n
e−1 S1 + e−1 un (−L) + 1 − e−1 un−1 (−L) −
u
(−L)
2
2

+e

−α(L−xj )

1 − 2e−1 n−2
5
˜n
u
(L) .
e−1 S2 + e−1 un (L) + 1 − e−1 un−1 (L) −
2
2

For the second order dissipative scheme, equation (5.29) now becomes

un+1 (x) =

L+ctn

α
−L−ctn

5 n
1
u (y) − un−1 (y) + un−2 (y) e−α|x−y| dy.
4
4

which in turn changes equations (5.30) to

un+1 (x) =
n+1 =
S1
n+1 =
S2

L

α
−L
−L

1
5 n
n+1
n+1
u (y) − un−1 (y) + un−2 (y) e−α|x−y| dy + S1 + S2 ,
4
4
α

1
5 n
u (y) − un−1 (y) + un−2 (y) e−α(x−y) dy,
4
4

α

5 n
1
u (y) − un−1 (y) + un−2 (y) e−α(y−x) dy.
4
4

−L−ctn
L+ctn
L

136

Repeating the same change of variables will give

˜n+1 = αc
S2

tn

5
1
u (L, tn − s) − u (L, tn−1 − s) + u (L, tn−2 − s) e−αcs ds
4
4

0

˜n
= e−αc∆t S2 + αc
√

˜n
= e− 2 S2 + αc

∆t
0
∆t

0

5
1
u (L, tn − s) − u (L, tn−1 − s) + u (L, tn−2 − s) e−αcs ds
4
4
1
5
u (L, tn − s) − u (L, tn−1 − s) + u (L, tn−2 − s) e−αcs ds.
4
4

Repeating the estimation of the linearized integral (5.32) now gives
√

√

2s
−
2
αc
u (L, tn
=
(zs + b)e ∆t ds
∆t 0
0 

√
√
√
2s ∆t
2s
∆t ∆t −
∆t −
2

−(zs + b) √ e ∆t
+
z √ e ∆t ds
=

∆t
2
2
0
0


√
√
2s ∆t
s ∆t
2
∆t −
∆t2 −

=
−
ze ∆t
−(zs + b) √ e ∆t

∆t
2
2
0
0
∆t

∆t

− s) e−αcs ds

√

∆t2 −√2
∆t
∆t2
∆t √
ze
+ b√ + z
−(z∆t + b) √ e− 2 −
2
2
2
2
√
√
∆t
∆t
= −(z∆t + b)e− 2 − √ ze− 2 + b + z √ .
2
2
=

2
∆t

Substituting back in for z and b from (5.31) produces
√

un−1 (L) − un (L) e− 2
un−1 (L) − un (L)
√
√
+ un (L) +
= −un−1 (L)e− 2 −
2
2
√
√
√
√ −√ 2
2 − 1 + e− 2 n
1 − 2e
− e− 2 n−1
√
√
=
u (L) +
u
(L)
2
2
√

= 1 un (L) + 2 un−1 (L)

137

with
√

√

2 − 1 + e− 2
√
=
and
1
2
√
√
√
1 − 2e− 2 − e− 2
√
.
2=
2

Thus
√
˜
˜n+1 = e− 2 S n + 5 1 un (L) +
S2
2
4

5
− 1 un−1 (L) −
4 2

2

+

1
n−2 (L) + 1 un−3 (L).
1 u
4
4 2

˜n+1 will be similar.
S1
Therefore, the fully discretized scheme will be given by
m

un+1
j

Aji

=
i=1

1
5 n
−α(L+xj ) ˜n+1
−α(L−xj ) ˜n+1
ui − un−1 + un−2 + e
S1 + e
S2 .
i
4
4 i

Similarly, for the second order non-dissipative scheme,

un+1 (x) + un−1 (x) = α
=α

L+ctn

e−α|x−y| un (y)dy

−L−ctn
L
e−α|x−y| un (y)dy
−L

n
n
+ S1 + S2 .

We then follow exactly the same procedure outlined for the second order dissipative
scheme to get
√

˜
˜n+1 = e− 2 S n + 1 un (L) + 1 un−1 (L)
S2
2

˜n+1
with a a similar result for S1 .
138

Therefore,
N

un+1
j

=

n−1
−uj

−α(L+xj ) ˜n
S1

Aji un + e
i

+

+e

−α(L−xj ) ˜n
S2 .

(5.32)

i=1

5.2.1

Performance

Instituting this non-reflecting boundary condition definitely produces an increase in performance as can be seen in table 5.4, based on the following parameters:
2

• The initial conditions are f (x) = e−100x and g(x) = 0.
• The domain is [−1, 1].
• ∆x = ∆t = 0.01 with a wave speed of c = 1.
• The 2-norm and ∞-norm were taken at time t = 2k, k ∈ {0, 1, 2, 3, 4, 5} as to measure
based on recombinations.
Implementing a non-reflecting boundary condition produces significantly better results
than an absorbing layer, especially when considering that the non-reflecting boundary needs
neither additional mesh points or an implicit solve. However, it does have certain drawbacks
that will force continued consideration of PML as a viable option for further work. First,
the non-reflecting boundary condition established here is strictly one dimensional. The idea
translates well to higher dimensions, but demands a considerably greater amount of work
and based on the literature, will produce lesser results. PML, on the other hand, was initially
created for a two dimensional setting and has been well established in three dimensions as
well.

139

Further, the non-reflecting boundary may not be stable under certain perturbations of the
problem. When first implementing this code, the mesh points xi were located at −L+i∗∆x,
meaning the domain truly covered [−L − ∆x/2, L + ∆x/2]. This actually caused the code to
diverge as exponential growth was seen at the boundary. When the domain was reduced to
[−L, L] and the boundary points were estimated by using the end mesh points ±(L − ∆x/2),
the code converged. Currently, a Taylor approximation is used for the value at the boundary
based upon the two mesh points nearest the boundary.

Table 5.4: Performance of the Dissipative 1st Order Scheme with a Non-reflecting Boundary
T ime
0
2
4
6
8
10

2-norm
3.54022
0.072160
0.0017082
0.000043195
0.0000011799
0.00000039215

140

∞-norm
1.00000
0.013608
0.00024755
0.00000047568
0.00000012592
0.0000000051711

Chapter 6
MOLT in Higher Dimensions Utilizing
an Alternating Directly Implicit Scheme
The adaptation of the scheme to higher dimensions creates a more difficult, yet more interesting challenge. Further, the MOLT scheme will be more valuable in higher dimensions
as the true known solution for the wave equation is considerably more complex in multiple
dimensions. Some of the difficulties that will arise due to higher dimensionality include an exponential growth in the number of mesh points as well as Green’s functions integrals that can
no longer can be evaluated analytically. Additionally, an increase in geometrical complexity
is to be expected as higher dimensions are more inclined to non-Cartesian coordinates, such
as spherical coordinates for three dimensions.
While there is currently no way to avoid the increase in the number of mesh points,
we can avoid the latter difficulties by implementing an Alternating Direct Implicit (ADI)
scheme, which will split the problem dimensionally and treat it as a series of one dimensional
problems.

141

6.1

Freespace Solution

Consider the wave equation (1.1) specifically in two dimensions with no applied boundary
conditions. That is,
∂ 2v
∂t2

= c2 2 v,

x ∈ R2 ,

v(x, 0) = f (x),

t ≥ 0,

vt (x, 0) = g(x).

Once again, the Method of Lines Transpose begins by discretizing the temporal derivative,
yielding
v n+1 − 2v n + v n−1
∆t2

= c2 2 v.

This can then be re-arranged into a Helmholtz type equation as

2

− α2 v =

−2v n + v n−1
c2 ∆t2

,

α=

1
.
c∆t

(6.1)

From here, note that this Helmholtz operator may be split dimensionally as

2

− α2 =
=
=

∂ ∂
∂ ∂
+
− α2
∂x ∂x ∂y ∂y
∂2

−1
α2

2
∂x

−1

∂2

α2

2
∂x

− α2
− α2

∂2
2
∂y

∂2
2
∂y

1 ∂2 ∂2

− α2

+

− α2

+ O ∆t2 .

2 2
α2 ∂x ∂y

While this dimensional split will introduce an O ∆t2 error if we drop the

142

(6.2)

1 ∂2 ∂2
2 2
α 2 ∂x ∂y

term,

this is acceptable as none of the schemes we are implementing are of order greater than 2.
While this may affect the accuracy of the scheme, it will not affect the order of the scheme.
We then take the split operator (6.2) and substitute it into equation (6.1) to obtain

−1

∂2

α2

2
∂x

∂2

− α2

2
∂y

− α2

v n+1

=

−2v n + v n−1
c2 ∆t2

.

This may then be arranged as the two step system of Helmholtz equations

Lx Ly v n+1 = −α2

−2v n + v n−1
c2 ∆t2

,

(6.3)

where

Lx =

∂2
2
∂x

− α2 and Ly =

∂2
2
∂y

− α2 .

Inverting the operators in order in equation (6.3) using our MOLT solver gives

Ly v n+1 (x, y) =
=⇒ v n+1 (x, y) =
= α4

−α2

−2v n (z) + v n−1 (z) −α|x−z|
e
dz
c2 ∆t2
−2v n (z) + v n−1 (z) −α|x−z|
−α2
e
dz e−α|y−z| dz
2 ∆t2
c
2v n (z) − v n−1 (z) e−α|x−z| dz e−α|y−z| dz.

This is the method for which we will implement the first order dissipative method. In

143

exactly the same fashion, we may calculate the second order dissipative scheme, given by

5 n
1
v (z) − v n−1 (z) + v n−2 (z) e−α|x−z| dz e−α|y−z| dz
4
4

v n+1 (x, y) = α4

and the second order non-dissipative scheme, given by

v n+1 (x, y) = −v n−1 (x, y) + α4

2v n (z)e−α|x−z| dz e−α|y−z| dz.

To minimize any error that may arise from always propagating first in one direction,
we will solve the problem twice, changing only the order of the dimensions. That is, find
numerical solutions to both

Lx Ly v n+1 = f v n , v n−1 , v n−2

and

Ly Lx v n+1 = f v n , v n−1 , v n−2

and take v n+1 (x, y) to be the average of these two solutions. While this is helpful for two
dimensions, solving for all possible orders of dimensions may become cost prohibitive as the
number of dimensions increases and n! solutions are needed for a domain in Rn .

6.2

Finite Domain

Thanks to the splitting of the operator in (6.2), enforcing any of the boundary conditions from
chapters 4 and 5 is as simple as applying the boundary condition in each step corresponding
to a dimension.
For example, to apply one-way wave boundary conditions to the 1st order dissipative
144

method, apply the temporal integral from 5.2 to each of the dimensions in equation (6.3),
which can be rewritten as

α2 2v n (z) − v n−1 (z) e−α|x−z| dz,

v=
˜

v n+1 = α2

v=
˜

∂2
2
∂y

− α2

v (z)e−α|y−z| dz.
˜

v n+1 ,

(6.4)
(6.5)

Thus, the fully discretized solution for (6.4) for each row u, upon which y is held constant,
˘
will be given by
m

˜
u=
˘
i=1

1
Aji un − un−1
˘i
˘
2 i

−α(L+xj )

1 − 2e−1 n−2
5
˜n
u
˘
(−L, y)
˘
e−1 Sx1 + e−1 un (−L, y) + 1 − e−1 un−1 (−L, y) −
˘
2
2

−α(L−xj )

5
1 − 2e−1 n−2
˜n
e−1 Sx2 + e−1 un (L, y) + 1 − e−1 un−1 (L, y) −
˘
˘
u
˘
(L, y) ,
2
2

+e
+e

˜
˜
˜
˜
where Sx1 and Sx2 are the natural equivalents to the original definition of S1 and S2 taken
at the boundaries x = −L and x = L, respectively.
This will then give that the fully discretized solution for (6.5) is
m

un+1 (x, y)

=

−α(L+yj ) −1 ˜n
˜
Aji ui + e
˘
e Sy1 + e−1 un (x, −L)

i=1
−1
5 −1
n−1 (x, −L) − 1 − 2e un−2 (x, −L)
+ 1− e
u
2
2
−α(L−yj )

+e

5
1 − 2e−1 n−2
˜n
e−1 Sy2 + e−1 un (x, L) + 1 − e−1 un−1 (x, L) −
u
(x, L) .
2
2

Any of the previous boundary conditions may be applied to any of the three schemes
already presented in a similar manner.
145

An example of the first order dissipative method conjoined with ADI with the one way
wave boundary condition of section 5.2 is shown in figures 6.1 - 6.3. The initial pulse was
−5 x2 +y 2

given by f (x, y) = e

and allowed to propagate outward at a speed of c = 1. The

wave was then tracked at times t = {0, 0.5, 1.25}.

Figure 6.1: Propagation of the First Order Dissipative ADI Scheme at t = 0

146

Figure 6.2: Propagation of the First Order Dissipative ADI Scheme at t = 0.5

Figure 6.3: Propagation of the First Order Dissipative ADI Scheme at t = 1.25

147

Chapter 7
Soft Sources
When computing the wave equation, a soft source is a prescribed source which is time varying
and introduces propagating waves into the medium at a given spatial location. The source
is said to be soft since, while it will emit waves at that given location, the source does not
interact with incoming waves and thus creates no reflection. To implement a soft source for
the one dimensional wave equation, begin with the equation (1.1) with zero initial conditions
which satisfies the soft source condition

v(xs , t) = S(t).

(7.1)

This soft source condition is not a boundary condition as the domain is free to extend to
the left and right of xs and waves may freely travel in either direction. Although we are
interested in time harmonic sources, this formulation will be for a general temporal function
S(t). We wish to understand how this problem might be reformulated by

vtt = c2 vxx + Q(x, t),

t > 0,
148

−∞ < x < ∞

(7.2)

with zero initial conditions, without explicitly enforcing condition (7.1). That is, we would
like to find a source term Q(x, t) that ensures v(xs , t) = S(t) will hold. We proceed by
constructing the Laplace solution to equation (1.1) with the soft source condition (7.1) and
equation (7.2) independently and equating the results. After applying the Laplace transform
to the first equation, it is easily seen that
s 2
V = Vxx ,
c

−∞ < x < ∞,

ˆ
and after imposing V (xs , s) = S(s), along with proper decay at ∞, this produces
s
ˆ
V (x, s) = S(s)e c

− |x−xs |

.

(7.3)

Likewise, after applying the Laplace transform to the second equation,
s 2
ˆ
V = Vxx + Q(x, s),
c

−∞ < x < ∞,

which can be solved using the Green’s function for the Helmholtz operator. Thus
s
∞
− |x−y|
c
ˆ
Q(y, s)e c
dy.
V (x, s) =
2s −∞
ˆ
We now impose upon Q(y, s) that equations (7.3) and (7.4) be equal, yielding
s
s
∞
− |x−y|
− |x−xs |
c
ˆ
ˆ
Q(y, s)e c
dy = S(s)e c
.
2s −∞

149

(7.4)

2s ˆ
ˆ
This may be solved by separating the variables in Q(y, s) as
S(s)g(y). This will require
c
that the undetermined function g(y) satisfies
∞

s
g(y)e c

− |x−y|

s
dy = e c

− |x−xs |

.

−∞

While the function g(y) may not be unique, there is the natural choice of

g(y) = δ(y − xs ),

which ensures that the solutions to equations (7.3) and (7.4) are equal. This gives that

2s ˆ
ˆ
Q(x, s) = S(s)δ(x − xs )
c

or, after taking the inverse Laplace transform,

2
Q(x, t) = S (t)δ(x − xs ).
c

Therefore, solving the wave equation (1.1) with the soft source condition (7.1) is equivalent
to solving the problem

2
utt = c2 uxx + S (t)δ(x − xs ),
c

150

t > 0, −∞ < x < ∞.

(7.5)

7.1

Implementation

The implementation of a soft source based on equation (7.5) follows in a manner similar to
the methods already presented. For the first order dissipative scheme, first discretize the
temporal derivative and then rearrange to form the Helmholtz equation

∂xx −

1
(c∆t)2

un+1

=

−2un (x) + un−1 (x)
(c∆t)2

+

2
c3

S (t)δ(x − xs ).

Upon inversion of the Helmholtz operator, we have

un+1 = α
=α

1
un (y) − un−1 (y) e−α|x−y| dy −
2

1
αc3

S (t)δ(y − xs )e−α|x−y| dy

(7.6)

1
1
un (y) − un−1 (y) e−α|x−y| dy −
S (t)e−α|x−xs | .
3
2
αc

Some care should be taken in choosing the source function S(t). Ideally, we would prefer
a source term that creates a wave and thus should be some combination of sin and cos
terms. Additionally, we would like the source to initially be off, i.e., S(0) = 0 and thus,
we will consider source terms of the from S(t) = sin (ω(t − (x − xs )). Further, we would
like to gradually introduce the source term in order to limit numerical error in the initial
introduction of the term; therefore, the source term will look like

S(t) = E(σt) sin (ω(t − (x − xs )),

where E(t) satisfies E(0) = 0, E(∞) = 1 and E ∈ C∞ . A good example of this is to
choose E(t) = tanh(t). The constants σ and ω are chosen to ensure that E(σt) and
sin (ω(t − (x − xs )) have different frequencies. Finally, if we wish to turn the signal off,
151

we can multiply S(t) by E2 (t), where E2 (t) satisfies E2 (0) = 1, E2 (∞) = 0 and E2 ∈ C∞ .
On a finite domain, we may include boundary conditions by simply adding in the boundary correction term, as was done in section 6.2. However, if we choose to utilize boundary
conditions, it makes more sense to include outflow boundary conditions as to be able to avoid
the continuous buildup of reflected waves that would occur with non-outflow boundary conditions.
The implementation of a soft source may be seen in figure 7.1. It simulates the first order
dissipative scheme with the one way wave boundary conditions from section 5.2. The wave
is created in the center of the domain using the source term S(t) = sin(πt) tanh(π/4). The
waves shown were taken at times t = {1, 3, 5} and show the general introduction of the wave.
The long term behavior, showing the waves flowing out of the domain is displayed in figure
7.2.
Wave Introduction via Soft Source
0.5

t=1
t=3
t=5

0.4
0.3
0.2
0.1
0
−0.1
−0.2
−0.3
−0.4
−5

−4

−3

−2

−1

0
x

1

2

Figure 7.1: Wave Introduction via Soft Source

152

3

4

5

Long Term Time Behavior of the Soft Source
0.6

t = 21
t = 32

0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−5

−4

−3

−2

−1

0
x

1

2

3

4

5

Figure 7.2: Long Term Time Behavior of the Soft Source

7.2

Higher Dimensions

When we moved the MOLT solver to a higher dimension, ADI was incorporated to split
the problem dimensionally. Thus, the natural choice for creating a soft source in a domain
encapsulating more than one dimension would be to split the source dimensionally as well.
The latter term from equation (7.6) will be applied to each dimension as it is split in the
ADI framework.
The simplest case of the implementation of a two-dimensional soft source in a dimensionally split environment is given as such. Begin with the free-space split problem, given
by equation (6.3) as

Lx Ly v n+1 = −α2

153

−2v n + v n−1
c2 ∆t2

.

In the first step, invert the operator in the x-dimension, Lx as well as propagating the wave
emanating from my source, as given by equation (7.6). This results in

w(x, y) =

Ω

α3

−2v n (z, y) + v n−1 (z, y) −α|x−z|
1
e
dz −
S (t)e−α|x−xs | ,
2
αc3

w(x, y) = Ly v n+1 (x, y).

Repeating the process in the y-dimension will then give that

v n+1 (x, y) =

Ω

−

1
w(x, z) −α|y−z|
e
dz −
S (t)e−α|y−ys | .
3
2
αc

As in section 6.1, it may be useful to solve the problem twice, varying the dimension
order and averaging the solution to minimize any dimensional order error.

154

Chapter 8
Future Directions

8.1

Further Problems

As the implementation of d’Alembert’s solution to the one-dimensional wave equation is
fairly straight-forward, more interesting problems lie in higher dimensions. The results that
have been proved here are useful as a basis for extending the work already done in one
dimension to higher dimensions. We have used an ADI scheme to dimensionally split the
problem and lift the numerical solution from one dimension to two dimensions. However,
this approach may not be practical for interesting problems in plasma physics, which will
likely require the six dimensions of space and velocity-space. The usefulness of the MOLT
method may be determined by how well it can be scaled for use in a higher dimensional
setting.
We would like to apply this work and determine its usefulness in solving other PDEs.

155

The most natural extension would be to the heat equation,

vt = c 2 v,

v(x, 0) = f (x).

While the Method of Lines Transpose idea has been used to numerically solve the heat
equation [51], the convergence of the method was never proven. The extension of the proofs
of consistency and stability should be fairly straight forward to extend from the wave equation
to the heat equation.
One of the factors impeding the practicality of the Method of Lines Transpose scheme is
its limitation to second order. The natural extensions to a third order dissipative scheme and
a third and fourth order purely dispersive scheme have also been investigated. Alas, these
methods were not feasible due to stability issues. For the third order dissipative scheme,
one of the eigenvalues of the matrix At was slightly larger than one, leading to a long-term
instability. Due to its symmetric temporal structure, the fourth order purely dispersive
scheme required that all the eigenvalues of At multiply to one. Numerical testing showed
two of the eigenvalues to be complex with norm equal to one, but unfortunately, the other
two were shown to be real and of different magnitude, causing the method to be completely
unstable. The issue with the third order purely dispersive scheme is slightly different. The
1
third order approximation to the temporal derivative vtt about the point tn− 2 is given by

vtt =

v n+1 − v n − v n−1 + v n−2
.
2

156

A centered third order equivalence for vxx is given by

vxx =

n−2
n−1
n
n+1
−vxx + 9vxx + 9vxx − vxx
.
16

Thus, rearranging the wave equation (1.1) into a Helmholtz type equation would give

n
n−1
L v n+1 − v n − v n−1 + v n−2 = 8vxx + 8vxx ,

L(v) =

∂xx +

8
(c∆t)2

v.

Due to the Helmholtz operator having a positive sign, the local nature of the inverted operator is destroyed as the corresponding Green’s function now becomes complex and oscillatory
instead of a negative exponential, yielding this as an unfeasible scheme.
Defect correction schemes [20, 23, 24, 31, 57] may provide a way to overcome these order
constraints. In a defect correction setting, the differential equation is originally solved with
a low order based scheme. An integral equation for the error is then determined, yielding a
correction term which is then applied to the solution, which in turn will raise the order of
the numerical solution. These methods would have to be adapted to fit this second order
partial differential equation, but may provide an avenue to obtaining higher order numerical
solutions for the wave equation using the MOLT as its base scheme.

8.2

Numerical Work

In this dissertation, we have discussed the advantages of implementing the MOLT scheme
for the wave equation. As we have the ability to take a large CFL, we expect to be able

157

to outperform standard explicit methods, however, a comparison on the accuracy and computational complexity of this method against standard finite difference methods and finite
element methods should be done.
Additionally, some of the proofs of the numerical properties could be more completely
solved. Namely, in section 3.2, we showed the dispersion relationship for the semi-discrete
methods. Attempts have been made to show the dispersion relation for the fully discrete
methods, but to this point, have been unsuccessful in showing that the second order dissipative scheme actually achieves second order. This is something that we hope to be able to
show in the future. Also, the stability bounds from section 4.2 for the finite domain schemes
can be tightened down. All of the numerical tests that we have run for the eigenvalues of the
matrix At = Ap + Ah , defined in section 4.2, have shown the eigenvalues to be positive and
within the predetermined bounds. The expectation is that the method is fully stable and
we would like to find a way to prove this, getting rid of the dependence on the Bauer-Fike
theorem to determine a region on which stability is guaranteed.

158

BIBLIOGRAPHY

159

Bibliography
[1] S. Abarbanel, D. Gottlieb, and J.S. Hesthaven. Well-posed perfectly matched layers for
advective acoustics. J. Comput. Phys., 154(2):266–283, September 1999.
[2] S. Abarbanel, D. Gottlieb, and J.S. Hesthaven. Long time behavior of the perfectly
matched layer equations in computational electromagnetics. J. Sci. Comput., 17(14):405–422, 2002.
[3] Saul Abarbanel and David Gottlieb. A mathematical analysis of the pml method. J.
Comput. Phys., 134:357–363, 1997.
[4] Saul Abarbanel and David Gottlieb. On the construction and analysis of absorbing
layers in cem. Appl. Num. Math., 27:331–340, 1998.
[5] Bradley Alpert, Leslie Greengard, and Thomas Hagstrom. An integral evolution formula
for the wave equation. J. Comput. Phys., 162:536–543, 2000.
[6] Bradley Alpert, Leslie Greengard, and Thomas Hagstrom. Nonreflecting boundary
conditions for the time-dependent wave equation. J. Comput. Phys., 180:270–296, 2002.
[7] Daniel Appelo, Thomas Hagstrom, and Gunilla Kreiss. Perfectly matched layers for
hyperbolic systems: general formulation, well-posedness and stability. J. Appl. Math.,
67(1):1–23, 2006.
[8] Uri M. Ascher, Robert M. M. Mattheij, and Robert D. Russell. Numerical solution of
boundary value problems for ordinary differential equations, volume 13 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1995. Corrected reprint of the 1988 original.
[9] Sergey Asvadurov, Vladimir Druskin, Murthy Guddati, and Leonid Knizhnerman. On
optimal finite-difference approximation of pml. J. Numer. Anal., 41(1):287–305, 2003.
[10] Peter Bartello and Stephen Thomas. The cost-effectiveness of semi-lagrangian advection.
Mon. Weather Rev., 124:2883–2924, 1996.
[11] F.L. Bauer and E.T. Fike. Norms and exclusion theorems. Numer. Math, 2:137–144,
1960.
[12] E. Becache, S. Fauqueux, and P. Joly. Stability of perfectly matched layers, group
velocities and anisotropic waves. J. Comput. Phys., 188:399–433, 2003.
160

[13] Eliane Becache and Patrick Joly. On the analysis of berenger’s perfectly matched layers
for maxwell’s equations. Mathematical Modelling and Numerical Analysis, 36(1):87–119,
2002.
[14] Eliane Becache, Peter Petropoulos, and Stephen Gedney. On the long-time behavior of
unsplit perfectly matched layers. IEEE Trans. Antennas and Propagation, 52(5):1335–
1342, 2004.
[15] Jean-Pierre Berenger. A perfectly matched layer for the absorption of electromagnetic
waves. J. Comput. Phys., 114:185–200, 1994.
[16] Jean-Pierre Berenger. Three-dimensional perfectly matched layer for the absorption of
electromagnetic waves. J. Comput. Phys., 127:363–379, 1996.
[17] Jean-Pierre Berenger. Improved pml for the fdtd solution of wave-structure interaction
problems. IEEE Trans. Antennas and Propagation, 45(3):466–473, March 1997.
[18] A. Bermudez, L. Hervella-Nieto, A. Prieto, and R. Rodriguez. An exact bounded perfectly matched layer for time-harmonic scattering problems. J. Sci. Comput., 30, 2007.
[19] L. Bonaventura. An introduction to semi - lagrangian methods for geophysical scale
flows. Lecture Notes, ERCOFTAC Leonhard Euler Lectures, SAM-ETH Zurich, 2004.
[20] Anne Bourlioux, Anita Layton, and Michael Minion. High-order multi-implicit spectral
deferred correction methods for problems of reactive flow. J. Comput. Phys., 189(2):651–
675, 2003.
[21] Roman Chapko and Rainer Kress. Rothe’s method for the heat equation and boundary
integral equations. J. Integral Equations and Appl., 9:47–69, 1997.
[22] H. Cheng, L. Greengard, and V. Rokhlin. A fast adaptive multipole algorithm in three
dimensions. J. Comput. Phys., 155(2):468–498, 1999.
[23] Andrew Christlieb, Benjamin Ong, and Jing-Mei Qiu. Comments on high order integrators embedded within integral deferred correction methods. Comm. Appl. Math and
Comp. Sci., 4(1):27–56, 2009.
[24] Andrew Christlieb, Benjamin Ong, and Jing-Mei Qiu. Integral deferred correction methods constructed with high order runge-kutta integrators. Math. Comp., 79(270):761–783,
2010.
[25] R. Coifman, V. Rokhlin, and S. Wandzura. The fast multipole method for the wave
equation: A pedestrian prescription. IEEE Trans. Antennas and Propagation, 35(3):7–
12, 1993.
[26] F. Collino and P.B. Monk. Optimizing the perfectly matched layer. Comput. Methods
Appl. Mech. Engrg., 164:157–171, 1998.

161

[27] Richard Courant, Eugene Isaacson, and Mina Rees. On the solution of nonlinear hyperbolic differential equations by finite differences. Communications on Pure and Applied
Mathematics, 5:243–255, 1952.
[28] Richard Courant and Peter Lax. On nonlinear partial differential equations with two
independent variables. Communications on Pure and Applied Mathematics, 2:255–273,
1949.
[29] Stephen Davis. A rotationally biased upwind difference scheme for the euler equations.
J. Comput. Phys., 56:65–92, 1984.
[30] J. Diaz and P. Joly. A time domain analysis of pml models in acoustics. Comput.
Methods Appl. Mech. Engrg., 195:3820–3852, 2006.
[31] Alok Dutt, Leslie Greengard, and Vladimir Rokhlin. Spectral deferred correction methods for ordinary differential equations. BIT, 40(2):241–266, 2000.
[32] Bjorn Engquist and Andrew Majda. Absorbing boundary conditions for the numerical
simulation of waves. Math. Comp., 31(139):629–651, 1977.
[33] Bjorn Engquist and Andrew Majda. Radiation boundary conditions for acoustic and
elastic wave calculations. Commun. Pure Appl. Math, 32:313–358, 1979.
[34] Maruizio Falcone and Roberto Ferretti. Convergence analysis for a class of high-order
semi-lagrangian advection schemes. SIAM J. Numer. Anal., 35(3):909–940, June 1998.
[35] M.G.G. Foreman. An accuracy analysis of boundary conditions for the forced shallow
water equations. J. Comput. Phys., 64:334–367, 1986.
[36] Zydrunas Gimbutas and Vladimir Rokhlin. A generalized fast multipole method for
nonoscillatory kernels. SIAM J. Sci. Comput., 24(3):796–817 (electronic), 2002.
[37] Dan Givoli. Non-reflecting boundary conditions. J. Comput. Phys., 94(1):1–29, 1991.
[38] S.K. Godunov. A difference scheme for numerical computation of discontinuous solutions
of equations of fluid dynamics. Math. Sbornik, 47:271–306, 1959.
[39] A. Greenbaum, L. Greengard, and G.B. McFadden. Laplace’s equation and the dirichletneumann map in multiply connected domains. J. Comput. Phys., 105:267–278, 1993.
[40] L. Greengard, M.C. Kropinski, and A. Mayo. Integral equation methods for stokes flow
and isotropic elasticity in the plane. J. Comput. Phys., 125(2):403–414, 1996.
[41] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comput.
Phys., 73(2):325–348, 1987.
[42] Marcus Grote and Imbo Sim. Efficient pml for the wave equation. Preprint, 2010.
[43] Laurence Halpern. Absorbing boundary conditions for the discretization schemes of the
one-dimensional wave equation. Math. Comp., 38(158):415–429, 1982.
162

[44] Laurence Halpern and Lloyd Trefethen. Wide-angle one-way wave equations. J. Acoust.
Soc. Amer., 84:1397–1404, 1988.
[45] Amiram Harten, Peter Lax, and Bram Van Leer. On upstream differencing and godunovtype schemes for hyperbolic conservation laws. SIAM Review, 25(1):34–61, 1983.
[46] U. Hornung. A parabolic-elliptic variational inequality. Manuscripta Math., 39(2):155–
172, 1982.
[47] Fang Hu. On absorbing boundary conditions for linearized euler equations by a perfectly
matched layer. J. Comput. Phys., 129:201–219, 1996.
[48] J. Jia and J. Huang. Krylov deferred correction accelerated method of lines transpose
for parabolic problems. J. Comput. Phys., 227(3):1739–1753, 2008.
[49] D.J Jones, J.C South Jr., and E.B Klunker. On the numerical solution of elliptic partial
differential equations by the method of lines. J. Comput. Phys., 9(3):496–527, 1972.
[50] J. Kacur.
Method of Rothe in evolution equations, volume 1192.
Berlin/Heidelberg, 1986.

Springer

[51] Mary Catherine Kropinski and Bryan Quaife. Fast integral equation methods applied
for rothe’s method applied to the isotropic heat equation. Comp. Math. Appl.
[52] Peter Lax. Nonlinear hyperbolic equations. Communications on Pure and Applied
Mathematics, 6:231–258, 1953.
[53] M. Lentine, J.T. Gretarsson, and R. Fedkiw. An unconditionally stable fully conservative semi-lagrangian method. J. Comput. Phys., 230:2857–2895, 2011.
[54] Randall J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002.
[55] Peijun Li, Hans Johnston, and Robert Krasny. A Cartesian treecode for screened
Coulomb interactions. J. Comput. Phys., 228(10):3858–3868, 2009.
[56] A. Mazzia and F. Mazzia. High-order transverse schemes for the numerical solution of
PDEs. J. Comput. Appl. Math., 82(1-2):299–311, 1997.
[57] Michael Minion. Semi-implicit spectral deferred correction methods for ordinary differential equations. Commun. Math Sci., 1(3):471–500, 2003.
[58] I.M. Navon, B. Neta, and M.Y. Hussaini. A perfectly matched layer formulation for
the linearized shallow water equation models: The split equation approach. Monthly
Weather Review, 132(6):1369–1378, 2004.
[59] Haim Nessyahu, Eitan Tadmor, and Tamir Tassa. The convergence rate of godunov
type schemes. SIAM J. Numer. Anal., 31(1):1–16, 1994.

163

[60] Satish Reddy and Lloyd Trefethen. Stability of the method of lines. Numer. Math,
62:235–267, 1992.
[61] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions.
J. Comput. Phys., 86:414–439, 1990.
[62] J.S. Sawyer. A semi-lagrangian method of solving the vorticity advection equation.
Tellus, 15(4):336–342, 1963.
[63] Andrew Selle, Ronald Fedkiw, ByungMoon Kim, Yingjie Liu, and Jarek Rossignac. An
unconditionally stable maccormack method. J. Sci. Comput., 35:350–371, 2008.
[64] A. Sommerfeld. Lectures on Theoretical Physics. Academic Press, New York, 1964.
[65] Eric Sonnendrucker, Jean Roche, Pierre Bertrand, and Alain Ghizzo. The semilagrangian method for the numerical resolution of the vlasov equation. J. Comput.
Phys., 149:201–220, 1999.
[66] A. Staniforth and J. Cote. Semi-lagrangian integration schemes for atmospheric models
- a review. Mon. Weather Rev., 119:2206–2223, 1991.
[67] John Strikwerda. Finite Difference Schemes and Partial Differential Equations, volume 2. Springer, 2004.
[68] Lloyd Trefethen and Laurence Halpern. Well-posedness of one-way wave equations and
absorbing boundary conditions. Math. Comp., 47(176):421–435, 1986.
[69] J. Verwer and J.M. Sanz-Serna. Convergence of method of lines approximations to
partial differential equations. Computing, 33(3-4):297–313, 1984.
[70] Wen-Chyuan Yueh and Sui Sun Cheng. Explicit eigenvalues and inverses of tri-diagonal
toeplitz matrices with four perturbed corners. ANZIAM, pages 361–387, 2008.
[71] A. Zafarullah. Application of the method of lines to parabolic partial differential equations with error estimates. J. ACM, 17:294–302, 1970.
[72] Li Zhao and Andreas Cangellaris. A general approach for the development of unsplitfield time-domain implementations of perfectly matched layers for fdtd grid truncation.
IEEE Microwave and Guided Wave Letters, 6(5):209–211, 1996.

164