Notes on Basic Finite Dierence Techniques for Time Dependent

PDEs

July 2003

Contents
1 Basic Finite Dierence Techniques for Time Dependent PDEs 1
1.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Types of IVPs (by example) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Wave and \Wave-Like" (\Hyperbolic"): The 1-d Wave Equation . . . . . . . . . . . . 2
1.2.2 Diusion (\Parabolic"): The 1-d Diusion Equation . . . . . . . . . . . . . . . . . . . 2
1.2.3 Schrodinger: The 1-d Schrodinger Equation . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Some Basic Concepts, Denitions and Techniques . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Residual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.2 Truncation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.3 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.4 Consistency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.5 Order of an FDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.6 Solution Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.7 Relation Between Truncation Error and Solution Error . . . . . . . . . . . . . . . . . 4
1.3.8 Deriving Finite Dierence Formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Sample Discretizations / FDAs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4.1 1-d Wave equation with xed (Dirichlet) boundary conditions . . . . . . . . . . . . . . 5
1.4.2 1-d Diusion equation with Dirichlet boundary conditions . . . . . . . . . . . . . . . . 7
1.5 The 1-D Wave Equation in More Detail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6.1 Heuristic Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6.2 Von-Neumann (Fourier) Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7 Dispersion and Dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.8 The Leap-Frog Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.9 Error Analysis and Convergence Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.9.1 Sample Analysis: The Advection Equation . . . . . . . . . . . . . . . . . . . . . . . . 18
1.10 Dispersion and Dissipation in FDAs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

1 Basic Finite Dierence Techniques for Time Dependent PDEs

There are several good reference texts for this material. Among my personal favorites are [1]-[4].

1.1 Preliminaries
We can divide time-dependent PDEs into two broad classes:
1. Initial-value Problems (Cauchy Problems), spatial domain has no boundaries (either innite or
\closed"|e.g. \periodic boundary conditions"

1
 2. Initial-Boundary-Value Problems, spatial domain nite, need to specify boundary conditions
Note: Even if physical problem is really of Type 1, nite computational resources ! nite spatial domain
! approximate as Type 2; will hereafter loosely refer to either type as an IVP.
Working Denition: Initial Value Problem
 State of physical system arbitrarily (usually) specied at some initial time t = t0 .
 Solution exists for t  t0 ; uniquely determined by equations of motion (EOM) and boundary conditions
(BCs).
Issues in Finite Dierence (FD) Approximation of IVPs
 Discretization (Derivation of FDA's)
 Solution of algebraic systems resulting from discretization
 Consistency
 Accuracy
 Stability
 Converegence
 Dispersion / Dissipation
 Treatment of Non-linearities
 Computational cost|expect O(N ) work (N  number of \grid points" (discrete events at which
approximate solution is computed)

1.2 Types of IVPs (by example)

In the following three examples, u is always a function of one space and one time variable, i.e. u  u(x; t).
Such a problem is often referred to as \1-d" by numericists, the time dimension being implicit in this
nomenclature. I will also use the subscript notation for partial dierentiation, e.g. ut  @t u.
1.2.1 Wave and \Wave-Like" (\Hyperbolic"): The 1-d Wave Equation
utt = c2 uxx c 2 R; (1)
u(x; 0) = u0(x)
ut (x; 0) = v0 (x)
1.2.2 Diusion (\Parabolic"): The 1-d Diusion Equation
ut = uxx  2 R;  > 0: (2)
u(x; 0) = u0 (x)
1.2.3 Schrodinger: The 1-d Schrodinger Equation
i t = h + V (x; t) 2C (3)
2m xx

(x; 0) = 0 (x)

Note: Although (x; t) is complex in this case, we can rewrite 3 as a system of 2 coupled scalar, real-valued
equations.

2
1.3 Some Basic Concepts, Denitions and Techniques
We will be considering the nite-dierence approximation (FDA) of PDEs, and as such, will generally we
interested in the continuum limit, where the mesh spacing, or grid spacing, usually denoted h, tends to 0.
Because any specic calculation must necessarily be performed at some specic, nite value of h, we will also
be (extremely!) interested in the way that our discrete solution varies as a function of h. In fact, we will
always view h as the basic \control" parameter of a typical FDA. Fundamentally, for sensibly constructed
FDAs, we expect the error in the approximation to go to 0, as h goes to 0.
Let
Lu = f (4)
denote a general dierential system. For simplicity and concreteness, you can think of u = u(x; t) as a single
function of one space variable and time, but the discussion in this section applies to cases in more independent
variables (u(x; y; t); u(x; y; z; t)


etc.), as well as multiple dependent variables (u = u = [u1 ; u2;


; un]).
In (4), L is some dierential operator (such as @tt @xx ) in our wave equation example), u is the unknown,
and f is some specied function (frequently called a source function) of the independent variables.
Here and in section 1.9 it will be convenient to adopt a notation where a superscript h on a symbol indicates
that it is discrete, or associated with the FDA, rather than the continuum. (Note, however, that for simplicity
of presentation, we will not adopt this notation in much of the development below). With this notation, we
will generically denote an FDA of (4) by
Lh uh = f h (5)
where u is the discrete solution, f is the specied function evaluated on the nite-dierence mesh, and Lh
h h

is the nite-dierence approximation of L.

1.3.1 Residual
Note that another way of writing our FDA is
L h uh f h = 0 (6)
It is often useful to view FDAs in this form for the following reason. First, we have a canonical view of
what it means to solve the FDA|\drive the left-hand side to 0". Furthermore, for iterative approaches to
the solution of the FDA (which are common, since it may be too expensive to solve the algebraic equations
directly), we are naturally lead to the concept of a residual. The residual is simply the level of \non-
satisfaction" of our FDA (and, indeed, of any algebraic expression). Specically, if u~h is some approximation
to the true solution of the FDA, uh, then the residual, rh , associated with u~h is just
rh  Lh u~h f h (7)
This leads to the view of a convergent, iterative process as being one which \drives the residual to 0".
1.3.2 Truncation Error
The truncation error,  h , of an FDA is dened by
 h  Lhu f h (8)
where u satises the continuum PDE (4). We note that the form of the truncation error can always be com-
puted (typically using Taylor series) from the nite dierence approximation and the dierential equations.
1.3.3 Convergence
Assuming that our FDA is characterized by a single discretization scale, h, we say that the approximation
converges i
uh ! u as h ! 0: (9)
Operationally (i.e. in practice), convergence is clearly our chief concern as numerical analysts, particularly
if there is reason to suspect that the solutions of our PDEs are good models for real phenomena. We note
that this is believed to be the case for many interesting problems in general relativistic astrophysics|the
two black hole problem being an excellent example.

3
1.3.4 Consistency
Assuming that the FDA with truncation error  h is characterized by a single discretization scale, h, we say
that the FDA is consistent if
 h ! 0 as h ! 0: (10)
Consistency is obviously a necessary condition for convergence.
1.3.5 Order of an FDA
Assuming that the FDA is characterized by a single discretization scale, h, we say that the FDA is p-th order
accurate or simply p-th order if
lim
h!0
 h = O(hp ) for some integer p (11)

1.3.6 Solution Error

The solution error, eh , associated with an FDA is dened by
eh  u u h (12)
1.3.7 Relation Between Truncation Error and Solution Error
Is is common to tacitly assume that
 h = O(hp ) ! eh = O(hp )
This assumption is often warranted, but it is extremely instructive to consider why it is warranted and to
investigate (following Richardson 1910 (!) [5]) in some detail the nature of the solution error. We will return
to this issue in more detail in section 1.9.
1.3.8 Deriving Finite Dierence Formulae
The essence of nite-dierence approximation of a PDE is the replacement of the continuum by a discrete
lattice of grid points, and the replacement of derivatives/dierential operators by nite-dierence expressions.
These nite-dierence expressions (nite-dierence quotients) approximate the derivatives of functions at grid
points, using the grid values themselves. All of the operators and expressions we need can easily be worked
out using Taylor series techniques. For example, let us consider the task of approximating the rst derivative
ux(x) of a function u(x), given a discrete set of values uj  u(jh) as shown in Figure 1. As it turns out,
x = x + j ∆x = x + j h
j 0 0
∆x

j−1 j j+1

Figure 1: A one-dimensional, uniform nite dierence mesh. Note that the spacing, 4x = h, between
adjacent mesh points is constant. In the text we tacitly assume that the origin, x0 , of our coordinate system
is x0 = 0.
given the three values u(xj h); u(xj ) and u(xj + h), which we will denote uj 1 ; uj , and uj+1 respectively,
we can compute an O(h2 ) approximation to ux (xj )  (ux )j as follows. Taylor expanding, we have
uj 1 = uj h(ux)j + 21 h2 (uxx)j 61 h3 (uxxx)j + 24 1 h4 (u ) + O(h5 )
xxxx j
uj = uj
uj+1 = uj + h(ux)j + 21 h2 (uxx)j + 61 h3 (uxxx)j + 24
1 h4 (u ) + O(h5 )
xxxx j

4
We now seek a linear combination of uj 1 ; uj , and uj+1 which yields (ux )j to O(h2 ) accuracy, i.e. we seek
c ; c0 and c+ such that
c uj 1 + c0 uj + c+ uj+1 = (ux)j + O(h2 )
This results in a system of three linear equations for uj 1 ; uj , and uj+1 :
c + c0 + c+ = 0
hc + hc+ = 1
1 h2 c + 1 h2 c = 0
2 2 +
which has the solution
c = 21h
c0 = 0
c+ = + 21h
Thus our O(h2 ) nite dierence approximation for the rst derivative is
u(x + h) u(x h) = u (x) + O(h2 ) (13)
2h x

Note that it may not be obvious to you a priori, that the truncation error of this approximation is O(h2 ),
since a naive consideration of the number of terms in the Taylor series expansion which can be eliminated
using 2 values (namely u(x + h) and u(x h)) suggests that the error might be O(h). The fact that the O(h)
term \drops out" is a consequence of the symmetry, or centering of the stencil, and is a common theme in
such FDAs (which, naturally enough, are called centred dierence approximations).
Using the same technique, we can easily generate the O(h2 ) expression for the second derivative, which uses
the same dierence stencil as the above approximation for the rst derivative.
u(x + h) 2u(x) + u(x h) = u (x) + O(h2 ) (14)
h2 xx

Exercise: Compute the precise form of the O(h2 ) terms in expressions (13) and (14).

1.4 Sample Discretizations / FDAs

1.4.1 1-d Wave equation with xed (Dirichlet) boundary conditions
utt = uxx (c = 1) 0  x  1; t  0 (15)
u(x; 0) = u0 (x)
ut (x; 0) = v0 (x)
u(0; t) = u(1; t) = 0

We now introduce a discrete domain (uniform grid) (xj ; tn )|part of which is shown in Figure 2.
tn  n 4t ; n = 0; 1; 2;



xj  (j 1) 4x ; j = 1; 2;


J
unj  u(n 4t ; (j 1) 4x )
4x = (J 1) 1
4t =  4x   \Courant number"

5
 t n+1
∆x
n
∆t

x n−1

j−1 j j+1

Figure 2: Portion of uniform nite-dierence mesh (grid) for 1-d time-dependent problem. Note that the
spacings in both the spatial and temporal directions are constant

Note: When solving wave equations using FDAs, we will typically keep  constant when we vary 4x .
Thus, our FDA will always be characterized by single discretization scale, h.
4x  h
4t  h
(Also note the Fortran-style indexing of the spatial grid index (j = 1; 2;


) and the C-style indexing of
the temporal one (n = 0; 1;


). This is a particular convention which I, as a predominantly Fortran
programmer, nd convenient.)
n+1

n−1

j−1 j j+1

Figure 3: Stencil (molecule/star) for \standard" O(h2 ) approximation of (15).

FDA: \standard O(h )" 2

Discretized Interior equation:

  1 4t 2 (u ) n + O( 4t 4 )
( 4t ) 2
unj +1 2unj + ujn 1
= (utt ) nj + 12 tttt j

= (utt ) j + O(h )
n 2

  1 4x 2 (u ) n + O( 4x 4 )
( 4x ) 2
unj+1 2unj + unj 1
= (uxx) nj + 12 xxxx j

= (uxx) j + O(h2 )
n

6
Putting these two together, we get the O(h2 ) approximation
unj +1 2unj + ujn 1
unj+1 2unj + unj
4t 2 = 4x 2
1
j = 2; 3;


; J 1 (16)
Note that a scheme such as (16) is often called a three level scheme since it couples three \time levels" of
data (i.e. unknowns at three distinct, disrete times tn 1 ; tn ; tn+1 .
Discretized Boundary conditions:
un1 +1 = unJ +1 = 0
Discretized Initial conditions:
We need to specify two \time levels" of data (eectively u(x; 0) and ut (x; 0)), i.e. we must specify
u0j ; j = 1; 2;


; J
u1j ; j = 1; 2;


; J
ensuring that the initial values are compatible with the boundary conditions.
Note that we can solve (16) explicitly for unj +1:
 
unj +1 = 2unj ujn 1 + 2 unj+1 2unj + ujn 1
(17)
Also note that (17) is actually a linear system for the unknowns unj +1 ; j = 1; 2;


; J ; in combination with
the discrete boundary conditions we can write
A un = b+1
(18)
where A is a diagonal J  J matrix and un+1 and b are vectors of length J . Such a dierence scheme for
an IVP is called an explicit scheme.
1.4.2 1-d Diusion equation with Dirichlet boundary conditions
ut = uxx ( = 1) 0  x  1; t  0 (19)
u(x; 0) = u0 (x)
u(0; t) = u( 1; t) = 0
We will use same discrete domain (grid) as for the 1-d wave equation.

FDA: Crank-Nicholson
This scheme illustrates a useful \rule of thumb": Keep dierence scheme \centred"
 centred in time, centred in space
 minimizes truncation error for given h
 tends to minimize instabilities
Discretization of time derivative:
  1 1 4t 2 (u ) n+ 21 + O( 4t 4)
4t 1
unj +1 unj = (ut ) nj + 2 + 24 ttt j (20)
1
= (ut ) n+ 2
j
+ O( 4t 2 )

7
 n+1

n
j−1 j j+1

n + 1/2
centre−point of scheme: x = x , t=t
j

Figure 4: Stencil (molecule/star) for O(h2 ) Crank-Nicholson approximation of (19).

O(h2 ) second-derivative operator:

 
Dxx unj  4x 2
unj+1 2unj + unj+1 (21)
1 4x 2 @
Dxx = @xx + 12 xxxx + O ( 4x ) (22)
4

(Forward) Time-averaging operator, t :

 
t unj  12 unj +1 + unj = ujn+ 2 + 81 4t 2 (utt ) nj + 2 + O( 4t 4 )
1 1
(23)
 
t = I + 18 4t 2 @tt + O( 4t 4) (24)
t=tn+1=2

where I is the identity operator. Assuming that 4t = O( 4x ) = O(h), is is easy to show (exercise) that
h i 1
t Dxx unj = (uxx) nj + 2 + O(h2 )
Putting the above results together, we are lead to the (O(h2 )) Crank-Nicholson approximation of (19):
unj +1 unj h ni
4t =  t
Dxx uj (25)
Written out in full, this is
2 3
unj +1 unj 1 unj+1
+1
2unj +1 + unj +11 unj+1 2unj + unj
= 4 + 1
5 j = 2; 3;


; J 1 (26)
4t 2 4x 2 4x 2

We can rewrite (26) in the form

a+ unj+1
+1
+ a0 unj +1 + a unj +11 = bj j = 2; 3;


; J 1 (27)
where
a+  21 4x 2
a0  4t 1 + 4x 2

a  12 4x 2
bj  4t 1
4x 2

un + 1 4x un 2
+ unj

j 2 j +1 1

8
which, along with the BCs (un1 +1 = unJ +1 = 0), is again a linear system of the form
A un = b
+1

for the \unknown vector" un+1 . This time, however, the matrix A, is not diagonal, and the scheme is called
implicit|i.e. the scheme couples unknowns at the advanced time level, t = tn+1 .
Note that A is a tridiagonal matrix: all elements Aij for which j 6= i + 1; i or i 1 vanish. The solution
of tridiagonal systems can be performed very eciently using special purpose routines (such as DGTSV in
LAPACK [6]): specically, the operation count for solution of (26) is O(J ).
Also note that we can immediately write down the analogous scheme for the Schrodinger equation (3):
h  hD n i + V (x )  n
n+1 n

i j
4t =
j
2m t xx j j t j
(28)
In this case we get a complex tridiagonal system, which can also be solved in O(J ) time, using, for example,
the LAPACK routine ZGTSV.

1.5 The 1-D Wave Equation in More Detail

Recall our \standard" O(h2 ) discretization:
 
unj +1 = 2unj ujn 1 + 2 unj+1 2unj + unj 1
; j = 2; 3;


; J 1
un1 +1 = unJ +1 = 0
As we have discussed, to initialize the scheme, we need to specify u0j and u1j , which is equivalent (in the limit
h ! 0) to specifying u(x; 0) and ut (x; 0).
Before proceeding to a discussion of a \proper initialization", let us brie y digress and consider the continuum
case, and, for the sake of presentation, assume that we are considering a true IVP on an unbounded domain;
i.e. we wish to solve
utt = uxx 1<x<1 ; t0 (29)
As is well known, the general solution of (29) is the superposition of an arbitrary left-moving prole (v =
c = 1), and an arbitrary right-moving prole (v = +c = +1); i.e.
u(x; t) = `(x + t) + r(x t) (30)
where (see Figure 5)
` : constant along \left-directed" characteristics
r : constant along \right-directed" characteristics

This observation provides us with an alternative way of specifying initial values, which is often quite conve-
nient in practice. Rather than specifying u(x; 0) and ut (x; 0) directly, we can specify initial left-moving and
right-moving parts of the solution, `(x) and r(x). Specically, we set
u(x; 0) = `(x) + r(x) (31)
d` (x) dr (x)
ut (x; 0) = `0 (x) r0 (x)  dx (32)
dx
Returning now to the solution of the nite-dierenced version of the wave equation, it is clear that given
the initial data (31{32), we can trivially initialize u0j with exact values, but that we can only approximately
initialize u1j . The question then arises: How accurately must we initialize the advanced values so as to ensure
second order (O(h2 )) accuracy of the dierence scheme?

9
 : "left−directed" characteristics, x + t = constant , l(x + t) = constant
: "right−directed" characteristics, x − t = constant , r(x − t) = constant

Figure 5: Characteristics of the wave equation: uxx = utt . Signals (disturbances) travel along the charac-
teristics (dashed and dotted lines.)

A brief, heuristic answer to this question (which can be more rigorously justied) is as follows. We have
4x = O(h), 4t = O(h) and the FDA is O(h2 ). Since the scheme is O(h2 ), we expect that
uexact(x; t) uFD(x; t) = O(h2 )
for arbitrary, xed, FINITE t. However, the number of time steps required to integrate to time t is
O( 4t 1) = O(h 1 ). Thus, the per-time-step error must be O(h2 )=O(h 1 ) = O(h3 ), and, therefore, we
require
(uFD ) 1j = (uexact) 1j + O(h3 )
We can readily accomplish this using (1) Taylor series and (2) the equation of motion to rewrite higher time
derivatives in terms of spatial derivatives:
u1j = u0j + 4t (ut) 0j + 21 4t 2 (utt ) 0j + O( 4t 3 ) (33)
= u0j + 4t (ut) + 21 4t 2 (uxx) 0j + O( 4t 3 ) (34)
which, using results from above, can be written as
u1j = (` + r) j + 4t (`0 r0 ) j + 21 4t 2 (`00 + r00 ) j (35)

1.6 Stability Analysis

One of the most frustrating|yet fascinating|features of FD solutions of time dependent problems, is that
the discrete solutions often \blow up"|e.g. oating-point over ows are generated at some point in the
evolution. Although \blow-ups" can sometimes be caused by legitimate (!) \bugs"|i.e. an incorrect
implementation|at other times it is simply the nature of the FD scheme which causes problems. We are
thus lead to consider the stability of solutions of dierence equations (as well as their dierential-equation
progenitors).
Let us again consider the 1-d wave equation (15) and let us now remark that this is a linear, non-dispersive
wave equation, a consequence of which is the fact that the \size" of the solution does not change with time:
ku(x; t)k  ku(x; 0)k ; (36)
where k
k is an suitable norm, such as the L2 norm:
Z 1 =1 2

ku(x; t)k  u(x; t) dx

2
: (37)
0

10
We will use the property captured by (36) as our working denition of stability. In particular, if you
believe (36) is true for the wave equation, then you believe the wave equation is stable.
Fundamentally, if our FDA approximation converges, then we expect the same behaviour for the dierence
solution:
kunjk  ku0j k : (38)
Now, we construct our FD solution by iterating in time, generating
u0j ; u1j ; u2j ; u3j ; u4j ;



in succession, using the FD equation
 
unj +1 = 2unj ujn 1 + 2 unj+1 2unj + unj 1
:
As it turns out, we are not guaranteed that (38) holds for all values of   4t = 4x . In fact, for certain 
(all  > 1, as we shall see), we have
kunjk  kuj k ; 0

and for those , kunk diverges from u, even (especially!) as h ! 0|that is, the dierence scheme is unstable.
In fact, for many wave problems (including all linear problems), given that a FD scheme is consistent (i.e.
so that ^ ! 0 as h ! 0), stability is the necessary and sucient condition for convergence (Lax's theorem).
1.6.1 Heuristic Stability Analysis
Let us write a general time-dependent FDA in the form
un = G[un ] ;
+1
(39)
where G is some update operator (linear in our example problem), and u is a column vector containing
sucient unknowns to writen the problem in rst-order-in-time form. For example, if we introduce a new,
auxiliary set of unknowns, vj , dened by
vnj = ujn 1 ;
then we can rewrite the dierenced-wave-equation (16) as
 
unj +1 = 2unj vnj + 2 unj+1 2unj + unj 1
; (40)
vnj +1 = unj ; (41)
so with
un = [un; vn ; un ; vn ;


unJ; vnJ ] ;
1 1 2 2

(for example), (40-41) is clearly of the form (39). Equation (39) provides us with a compact way of describing
the solution of the FDA. Given initial data, u0 , the solution after n time-steps is
un = Gn u ; 0
(42)
where Gn is the n-th power of the matrix G. Now, assume that G has a complete set of orthonormal
eigenvectors
ek ; k = 1; 2;


J ;
and corresponding eigenvalues
k ; k = 1; 2;


J ;
so that
G ek = k ek ; k = 1; 2;


J :
11
We can then write the initial data as (spectral decomposition):
X
J
u =
0
c0k ek ;
k=1

where the c0k are coecients. Using (42), the solution at time-step n is then
X
J
!
u = G
n n
ck ek
0
(43)
k=1
X
J
= c0k (k )n ek : (44)
k=1

Clearly, if the dierence scheme is to be stable, we must have

jk j  1 k = 1; 2;


J (45)
p
(Note: k will be complex in general, so jj denotes complex modulus, jj  ? ).
Geometrically, then, the eigenvalues of the update matrix must lie on or within the unit circle (see Figure 6).
Im

unit circle

Re

Figure 6: Schematic illustration of location in complex plane of eigenvalues of update matrix G. In this
case, all eigenvalues (dots) lie on or within the unit circle, indicating that the corresponding nite dierence
scheme is stable.

1.6.2 Von-Neumann (Fourier) Stability Analysis

Von-Neumann stability analysis is based on the ideas sketched above, but additionally assumes that the
dierence equation is linear with constant coecients, and that the boundary conditions are periodic. We
can then use Fourier analysis, which has the same benets in the discrete domain|dierence operators
in real-space variable x ! algebraic operations in Fourier-space variable k|as it does in the continuum
Schematically, instead of writing
un+1 (x) = G[un(x)] ;

12
we consider the Fourier-domain equivalent:
u~ n (k) = G~ [~un (k)] ;
+1

where k is the wave-number (Fourier-space variable) and u~ and G ~ are the Fourier-transforms of u and G,
respectively. Specically, we dene the Fourier-transformed grid function via
Z 1
u~ n(k) = p12
+
e ikx
un (x) dx : (46)
1
For a general dierence scheme, we will nd that
u~ n (k) = G~ () u~ n(k) ;
+1

where   kh, and we will have to show that G ~ ( )'s eigenvalues lie within or on the unit circle for all
conceivable  . The appropriate range for  is
   ;
since the shortest wavelength representable on a uniform mesh with spacing h is  = 2h (Nyquist limit),
corresponding to a maximum wave number k = (2)= = =h.
Let us consider the application of the Von-Neumann stability analysis to our current model problem. We
rst dene a (non-divided) dierence operator D2 as follows:
D2 u(x) = u(x + h) 2u(x) + u(x h) :
Then, suppressing the spatial grid index, we can write the rst-order form of the dierence equation (40-41)
as
un+1 = 2un vn + 2 D2 un ;
vn+1 = un ;
or  u n +1  2+ D 1 2
u n: 2
 
v = 1 0 v (47)
In order to perform the Fourier transform, we need to know the action of D2 in Fourier-space. Using the
transform inverse to (46) we have
Z 1
u(x) = p1
+
eikx u~(k) dk ;
2 1
so
Z + 1
eikx u~(k) dk
D2 u(x) = u(x + h) 2u(x) + u(x h) = eikh 2 + e ikh
1
Z +1
eikx u~(k) dk :
= ei 2 + e i
1
Now consider the quantity 4 sin2 (=2):
 
4 sin2 2 = 4 e 2ie
i=2 i=2 2

 i=  2
= e 2
e i=2
= ei 2 + e i
;

13
so Z 1 
D2 u(x) = p1 4 sin2 2 eikx u~(k) dk :
+

2 1
In summary, under Fourier transformation, we have
u(x) ! u~ (k) ;
D2 u(x) ! 4 sin2 2 u~ (k) :
Using this result in the Fourier transform of (47), we see that we need to compute the eigenvalues of
 2 4 sin2 (=2) 1 ;
2

1 0
and determine the conditions under which the eigenvalues lie on or within the unit circle. The characteristic
equation (whose roots are the eigenvalues) is
 
 2 4 sin (=2)  1 
2 2
= 0
 1 ; 
or  

2 + 42 sin2 2 2  + 1 = 0 :
This equation has roots
    =
( ) = 1 22 sin2 2  1 22 sin2 2
1 2
1 :
We now need to nd sucient conditions for
j( )j  1;
or equivalently
j( )j  1:
2

To this end, we note that we can write

( ) = (1 Q)  ((1 Q)2 1)1=2 ;
where the quantity, Q
Q  2 sin2 2 ;
is real and non-negative (Q  0). There are now two cases to consider:
1. (1 Q)2 1  0 ,
2. (1 Q)2 1 > 0 .
In the rst case, ((1 Q)2 1)1=2 is purely imaginary, so we have
j( )j2 = (1 Q)2 + (1 (1 Q)2 ) = 1 :
In the second case, (1 Q)2 1 > 0 ! (1 Q)2 > 1 ! Q > 2, and then we have
1 Q ((1 Q2 ) 1)1=2 < 1 ;
so, in this case, our stability criterion will always be violated. We thus conclude that a necessary condition
for Von-Neumann stability is
(1 Q)2 1  0 ! (1 Q)2  1 ! Q  2 :

14
Since Q  2 sin2 (=2) and sin2 (=2)  1, we must therefore have
4t  1 ;
 4 x
for stability of our scheme (16). This condition is often called the CFL condition|after Courant, Friedrichs
and Lewy who derived it in 1928 (the ratio  = 4x = 4t is also frequently called the Courant num-
ber). In practical terms, we must limit the time-discretization scale , 4t , to values no larger than the
space-discretization scale, 4x . Furthermore, this type of instability has a \physical" interpretation, often
summarized by the statement the numerical domain of dependence of an explicit dierence scheme must
contain the physical domain of dependence.

1.7 Dispersion and Dissipation

Let us now consider an even simpler model \wave equation" than (1), the so-called advection, or color
equation:
ut = a ux (a > 0) 1<x<1 ; t0 (48)
u(x; 0) = u0 (x)
which has the exact solution
u(x; t) = u0 (x + at) (49)
Equation (48) is another example of a non-disspative, non-dispersive partial dierential equation.
To remind ourselves what \non-dispersive" means, let us note that (48) admits \normal mode" solutions:
u(x; t)  eik(x+at)  ei(kx+!t)
where !  ka; in general, of course, !  !(k) is known as the dispersion relation, and
d!  speed of propagation of mode with wave number k
dk
In the current case, we have
d! = a = constant
dk
which means that all modes propagate at the same speed, which is precisely what is meant by \non-
dispersive". Further, if we consider resolving the general initial prole, u0(x) into \normal-mode" (Fourier)
components, we nd that the magnitudes of the components are preserved in time, i.e. equation (48) is also
non-dissipative.
Ostensibly, we would like our nite-dierence solutions to have the same properties|i.e. to be dissipationless
and dispersionless, but, in general, this will not be (completely) possible. We will return to the issue of
dissipation and dispersion in FDAs of wave problems below.

1.8 The Leap-Frog Scheme

First note that (48) is a good prototype for the general hyperbolic system:
ut = Aux (50)
where u(x,t) is the n-component solution vector:
u(x; t) = [u (x; t); u (x; t);


un (x; t)]
1 2 (51)
and the n  n matrix A has distinct real eigenvalues
1 ; 2 ;


n

15
so that, for example, there exists a similarity transformation S such that
SAS = diag( ;  ;


n)
1
1 2

The leap-frog scheme is a commonly used nite-dierence approximation for hyperbolic systems. In the
context of our simple scalar (n = 1) advection problem (48):
ut = a ux
an appropriate stencil is shown in Figure 7. Applying the usual O(h2 ) approximations to @x and @t , our
n+1

n−1

j−1 j j+1

Figure 7: Stencil (molecule/star) for leap-frog scheme as applied to (48). Note that the central grid point
has been lled in this gure to emphasize that the corresponding unknown, unj , does not appear in the local
discrete equation at that grid point (hence the term \leap-frog")
leap-frog (LF) scheme is
unj +1 ujn 1
unj+1 unj
2 4t =a 2 4x
1
(52)
or explicitly  
unj +1 = ujn 1 + a unj+1 unj 1
(53)
where
4t
 4 x
is the Courant number as previously.
Exercise: Perform a von Neumann stability analysis of (52) thus showing that a  1 (which, you should
note, is just the CFL condition) is necessary for stability.
Observe that the LF scheme (52) is a three-level scheme. As in our treatment of the wave equation, utt = uxx
using the \standard scheme", we need to specify
u0j ; u1j j = 1; 2;


J
to \get the scheme going"|that is, we need to specify two numbers per spatial grid point. This should
be contrasted to the continuum where we need to specify only one number per xj , namely u0 (xj ). Again,
the initialization of the u0j is trivial, given the (continuum) initial data u0 (x), and, again, we need u1j to
O( 4t 3) = O(h3 ) accuracy for O(h2 ) global accuracy. Two possible approaches are as follows.
Taylor Series: The development here is parallel to that for the wave equation. We have
u1j = u0j + 4t (ut) 0j + 21 4t 2 (utt) 0j + O( 4t 2 )

16
also, from the equation of motion ut = aux, we get
utt = (ut )t = (aux)t = a (ut )x = a2 uxx:
so we have our desired initialization formula:


u1j = u0j + 4t (u00 ) 0j + 21 4t 2 a2 u000 0j + O( 4t 3) (54)

Self-Consistent Iterative Approach: The idea here is to initialize the u1j from the u0j and a version of the
discrete equations of motion which introduces a \cticious" half-time-level|see Figure 8.
1

1/2

0
j−1 j j+1

Figure 8: Stencil for initialization of leap-frog scheme for to (48). Note the introduction of the \ctitious"
half-time level t = t1=2 (squares).
Applying the leap-frog scheme on the stencil in the Figure, we have
1 1
u1j u0j uj2+1 uj2
4t = a 2 4x
1

or, explicitly solving for u1j :

 
u1j = u0j + 12  uj2+1 uj2
1 1
1

It is a straightforward exercise to show that in order to retain O(h2 ) accuracy of the dierence scheme,
we need \ctitious-time" values, u1j =2 which are accurate to O(h2 ) (i.e. we can neglect terms which are of
O(h2 )). In particular, if we dene u1j =2, via
1 u1j + u0j
u2
2
j
=
which amounts to dening the half-time values via linear interpolation in the advanced and retarded un-
knowns, we will retain second-order accuracy.
We are thus led to the following initialization algorithm which is perhaps best expressed in pseudo-code
(note, all loops over j are implicit:)
u[0,j] := u_0(x_j)
u[1,j] := u_0(x_j)
DO
usave[j] := u[1,j]
u[1/2,j] := (u[1,j] + u[0,j]) / 2

u[1,j] := u[0,j] + (lambda / 2) * (u[1/2,j+1] - u[1/2,j-1])

UNTIL norm(usave[j] - u[1,j]) < epsilon

17
1.9 Error Analysis and Convergence Tests
As a side remark, we note that the discussion in this section applies to essentially any continuum problem
which is solved using FDAs on a uniform mesh structure. In particular, the discussion applies to the
treatment of ODEs and elliptic problems, where, in fact, convergence is often easier to achieve due to the
fact that the FDAs are typically intrinsically stable (i.e. we have an easier time constructing stable FDAs for
these types of problems). We also note that departures from non-uniformity in the mesh do not, in general,
complete destroy the picture, but, rather, tend to distort it in ways which are beyond the scope of these
notes. In any case, my colleagues have been known to paraphrase my whole approach to this subject as
Convergence!, Convergence!, Convergence!
1.9.1 Sample Analysis: The Advection Equation
We again consider the solution of the advection equation, but this time we impose periodic boundary con-
ditions on our spatial domain, which we take to be 0  x  1 with x = 0 and x = 1 identied (i.e. we solve
the wave equation on S1  R):
ut = a ux (a > 0) 0  x  1; t  0 (55)
u(x; 0) = u0(x)
Note that the initial conditions u0 (x) must be compatible with periodicity, i.e. we must specify periodic
initial data.
Again, given the initial data, u0(x), we can immediately write down the full solution
u(x; t) = u0 (x + a t mod 1) (56)
where mod is the usual modulus function which \wraps" x + a t, t > 0 onto the unit circle. As we shall
see, because of the simplicity and solubility of this problem, one can perform a rather complete \analytic"
(i.e. closed form) treatment of the convergence of simple FDAs of (55). The point of the exercise, however,
is not to advocate parallel \analytic" treatments for more complicated problems. Rather, the key idea to
be extracted from the following is that, in principle (always), and in practice (almost always, i.e. I've never
seen a case where it didn't work, but then there's a lot of computations I haven't seen):
The error, eh, of an FDA is no less computable than the solution, uh itself.
This has widespread ramications, one of which is that there's really no excuse for publishing solutions of
FDAs without error bars, or their equivalents!
Proceeding with our sample error analysis of the leap-frog scheme applied to the advection equation, we rst
introduce some dierence operators for the usual O(h2 ) centred approximations of @x and @t :
unj+1 unj 1
Dx u n
j
 2 4x (57)
ujn+1
ujn 1
Dt unj  2 4t (58)
We again take
4x  h 4t   4x =  h
and will hold  xed as h varies, so that, as usual, our FDA is characterized by the single scale parameter,
h.
The idea behind our error analysis is that we want to view the solution of the FDA as a continuum problem,
and hence we will express both the dierence operators and the FDA solution as asymptotic series (in h) of

18
dierential operators, and continuum functions, respectively. We have the following expansions for Dx and
Dt :
Dx = @x + 16 h2 @xxx + O(h4 ) (59)
Dt = @t + 61 2 h2 @ttt + O(h4 ) (60)

Now, in terms of the general, abstract formulation of (1.3), we have:

Lu f = 0 () (@t a @x ) u = 0 (61)
L u fh = 0
h h
() (Dt a Dx) uh = 0 (62)
Lh u f h   h ()


(Dt a Dx) u   h = 61 h2 2 @ttt a @xxx u + O(h4 ) = O(h2 ) (63)
The Richardson ansatz: The key to our analysis is L.F. Richardson's old observation (ansatz) [5], that
the solution, uh , of any FDA which (1) uses a uniform mesh structure with scale parameter h, and (2) is
completely centred, should have the following expansion in the limit h ! 0:
uh(x; t) = u(x; t) + h2 e2 (x; t) + h4 e4 (x; t) +


(64)
Here u is the continuum solution, while e2 , e4 ,


are (continuum) error functions which do not depend on
h. In a very real sense (64), is the key expression from which all error analysis of FDAs derives. We note
that in the case that the FDA is not completely centred, we will have to modify the ansatz. In particular,
for rst order schemes (which are more common in relativistic astrophysics than one might expect!), we will
have
uh(x; t) = u(x; t) + he1(x; t) + h2 ex(x; t) + h3 e3 (x; t) +


(65)
Note that the Richardson ansatz (64) is completely compatible with the assertion discussed in (1.3.7), namely
that
 h = O(h2 ) ! eh  u uh = O(h2 ) (66)
However, the Richardson form (64) contains much more information than \second-order truncation error
should imply second-order solution error", telling us the precise form of the h dependence of uh.
Given the Richardson expansion, we can now proceed with our error analysis. We start from the FDA,
Lhuh f h = 0, and replace both Lh and uh with continuum expansions:
!


Lh uh = 0
(Dt a1Dx) u + h e2 +

1
= 0
2


! @t + 6 2 h2 @ttt a @x 6 ah2 @xxx +


u + h2 e2 +


= 0 (67)

We now demand that terms in (67) vanish order-by-order in h. At O(1) (zeroth-order), we have
(@t a @x) u = 0 (68)
which is simply a statement of the consistency of the dierence approximation. More interestingly, at O(h2 )
(second-order), we nd
(@t a @x ) e2 = 61 a@xxx 2 @ttt u

(69)
which, assuming that we view u as a \known" function, is simply a PDE for the leading order error function,
e2 . Moreover, the PDE governing e2 is of precisely the same nature as the original PDE (48).
In fact, we can solve (69) for e2 . Given the \natural" initial conditions
e2 (x; 0) = 0

19
(i.e. we initialize the FDA with the exact solution so that uh = u at t = 0), and dening q(x + at):


q(x + at)  61 a 1 2 a2 @xxx u(x; t)
we have
e2 (x; t) = t q(x + at mod 1) (70)
We note that, as is typical for leap-frog, we have linear growth of the nite dierence error with time (to
leading order in h). We also note that we can obviously push this analysis to higher order in h|what results,
then, is an entire hierarchy of dierential equations for u and the error functions e2, e4 , e6 ,


. Indeed, it is
useful to keep this view in mind:
When one solves an FDA of a PDE, one is not solving some system which is \simplied" realtive
to the PDE, rather, one is solving a much richer system consisting of an (innite) hierarchy of
PDEs, one for each function appearing in the Richardson expansion (64).
In the general case, of course, we will not be able to solve the PDE governing u, let alone that governing
e2 |otherwise we wouldn't be considering the FDA in the rst place! But it is precisely in this instance
where the true power of Richardson's observation is evident. The key observation is that starting from (64),
and computing FD solutions using the same initial data, but with diering values of h, we can learn a great
deal about the error in our FD approximations. The whole game of investigating the manner in which a
particular FDA converges or doesn't (i.e. looking at what happens as one varies h) is known as convergence
testing. It is important to realize that there are no hard and fast rules for convergence testing; rather, one
tends to tailor the tests to the specics of the problem at hand, and, being largely an empirical approach,
one gains experience and intuition as one works through more and more problems. That said, I emphasize
again that the Richardson expansion, in some form or other, always underlies convergence analysis of FDAs.
A simple example of a convergence test, and the one I use most often in practice is the following. We compute
three distinct FD solutions uh, u2h , u4h at resolutions h, 2h and 4h respectively, but using the same initial
data (as naturally expressed on the 3 distinct FD meshes). We also assume that the nite dierence meshes
\line up", i.e. that the 4h grid points are a subset of the 2h points which are a subset of the h points, so
that, in particular, the 4h points constitute a common set of events (xj ; tn ) at which specic grid function
values can be directly (i.e. no interpolation required) and meaningfully compared to one another. From the
Richardson ansatz (64), we expect:
uh = u + h2 e2 + h4 e4 +



u2h = u + (2h)2 e2 + (2h)4 e4 +



u4h = u + (4h)2 e2 + (4h)4 e4 +



We then compute a quantity Q(t), which I will call a convergence factor, as follows:

Q(t)  kkuu2h uuhkkx

4h 2h
(71)
x

where k
kx is any suitable discrete spatial norm, such as the `2 norm, k
k2 :

0 J 1=
X

1 2

kuhk = @J
2 uhj A
1 2
(72)
j =1

and, for concreteness, the subtractions in (71) can be taken to involve the sets of mesh points which are
common between u4h and u2h, and between u2h and uh . It is a simple exercise to show that, if our nite
dierence scheme is converging, then we should nd:
lim
h!0
Q(t) = 4: (73)

20
In practice, one can use additional levels of discretization, 8h, 16h, etc. to extend this test to look for \trends"
in Q(t) and, in short, to convince oneself (and, with luck, others), that the FDA really is converging.
Moreover, once convergence of an FDA has been established, then a point-wise subtraction of any two
solutions computed at dierent resolutions, will immediately provide an estimate of the level of error in
both. For example, if we have uh and u2h, then, again by the Richardson ansatz we have




u2h uh = u + (2h)2 e2 +


u + h2 e2 +


= 3h2e2 + O(h4 )  3eh  34 e2h (74)

Richardson extrapolation: Richardson's observation (64) also provides the basis for all the techniques of
Richardson extrapolation, where solutions computed at dierent resolutions are linearly combined so as to
eliminate leading order error terms, and hence provide more accurate solutions. As an example, given uh
and u2h which satisfy (64), we can take the linear combination, uh:
uh  4u 3 u
h 2h
(75)
which, by (64), is easily seen to be O(h4 ), i.e. fourth-order accurate!
4 u + h2 e2 + h4 e4 +




u + 4h2e + 16h4e +




u 
h 4u h u 2h
= 2 4

3 3
= 4h4e4 + O(h6 ) = O(h4 ) (76)
When it works, Richardson extrapolation has an almost magical quality about it, but one generally has to
start with fairly accurate (on the order of a few %) solutions in order to see the dramatic improvement in
accuracy suggested by (76). Partly because it is still a struggle to achieve that sort of accuracy (i.e. a few
%) for any computation in many areas of numerical relativity/astrophysics, techniques based on Richardson
extrapolation have not had a major impact in this context.
Independent Residual Evaluation A question which often arises in discussions of convergence testing is the
following:
\OK, you've established that uh is converging as h ! 0, but how do you know you're converging
to u, the solution of the continuum problem?"
Here, the notion of an independent residual evaluation is very useful. The idea is as follows: we have our
continuum PDE
Lu f = 0 (77)
and our FDA
L h uh f h = 0 (78)
We have demonstrated that u is apparently converging by, for example, computing the convergence fac-
h

tor (71) and verifying that it tends to 4 as h tends to 0. However, we do not know if we have derived and/or
implemented our discrete operator Lh correctly. Note that implicit in the \implementation" is the fact that,
particularly for multi-dimensional and/or implicit and/or multi-component FDAs, considerable \work" (i.e.
analysis and coding) may be involved in setting up and solving the algebraic equations for uh . As a check
that we are converging to u, we consider a distinct (i.e. independent) discretization of the PDE:
L~ hu~h f h = 0 (79)
The only thing we need from this FDA for the purposes of the independent residual test is the new FD
operator L~ h. As with Lh, we can expand L~ h in powers of the mesh spacing:
L~ h = L + h2 E2 + h4 E4 +


(80)
where E2 , E4 ,


are higher order (involve higher order derivatives than L) dierential operators. We then
simply apply the new operator L~ h to our FDA uh and investigate what happens as h ! 0. If uh is converging
to the continuum solution, u, we will have
uh = u + h2 e2 + O(h4 ) (81)

21
and we will compute



L~ h uh = L + h2 E2 + O(h4 ) u + h2 e2 + O(h4 ) = Lu + h2 (E2 u + L e2) = O(h2 ) (82)
i.e., L~ huh will be a residual-like quantity which converges quadratically as h ! 0. Conversely, if we have
goofed in our derivation and/or implementation of Lhuh = f h = 0, but we still see convergence; i.e. we
have, for example, u2h uh ! 0 as h ! 0, then we must have something like
uh = u + e0 + he1 + h2 e2 +


(83)
where the crucial fact is that the error must have an O(1) component, e0 . In this case, we will compute



L~ h uh = L + h2 E2 + O(h4 ) u + e0 + he1 + h2 e2 + O(h4 ) = Lu + Le0 + hLe1 + O(h2 ) = Le0 + O(h) (84)
and, unless we are extraordinarily lucky, and L e0 vanishes, we will not observe the expected convergence,
rather, we will see L~ huh f h tending to a nite (O(1)) value|a sure sign that something is wrong.
There is of course, the problem that we might have slipped up in our implementation of the \independent
residual evaluator", L~ h , in which case the results from our test will be ambigous at best! However, a key
point here is that because L~ h is only used a posteriori on a computed solution (we never use it to compute
u~h, for example) it is a relatively easy matter to ensure that L~ h has been implemented in an error-free fashion
(perhaps using symbolic manipulation facilities). Furthermore, many of the restrictions commonly placed
on the \real" discretization (such as stability and the ease of solution of the resulting algebraic equations)
do not apply to L~ h .

1.10 Dispersion and Dissipation in FDAs

We again consider the advection model problem, ut = a ux, but now discretize only in space (semi-
discretization) using the usual O(h2 ) centred dierence approximation:
uj+1 uj
ut = a Dx u  a 2 4x
1
(85)
We now look for normal-mode solutions to (85) of the form
u = eik (x+a t)
0

where the \discrete phase speed", a0 , is to be determined. Substitution of this ansatz in (85) yields
ika0 u = a (2i sin(k 4x )) u
2 4x
or, solving for the discrete phase speed, a0
a0 = a sin(k 4x ) = a sin 
k 4x 
where we have dened the dimensionless wave number,  :
  k 4x
In the low frequency limit,  ! 0, we have the expected result:
a0 = a sin  ! a

so that low frequency components propagate with the correct phase speed, a. However, in the high frequency
limit,  ! , we have
a0 = a sin  ! 0 !!

22
i.e. the highest frequency components of the solution don't propagate at all! This is typical of FDAs of wave
equations, particularly for relatively low-order schemes. The propagation of high frequency components of
the dierence solution is essentially completely wrong. Arguably then, there can be little harm in attenuating
(dissipating) these components, and, in fact, since high frequency components are potentially troublesome
(particularly vis a vis non-linearities and the treatment of boundaries), it is often advantageous to use a
dissipative dierence scheme.
Some FDAs are naturally dissipative (the Lax-Wendro scheme, for example), while others, such as leap-
frog, are not. In the case of a leap-frog-based scheme, the idea is to add dissipative terms to the method,
but in such a way as to retain O(h2 ) accuracy of the scheme. Consider, for example, the leap-frog scheme
as applied to the advection model problem:
 
unj +1 = ujn 1 + a unj+1 unj 1

We add dissipation to the scheme by modifying it as follows:

  n 
unj +1 = ujn 1 + a unj+1 unj 1 16 uj +2
1
4ujn+11 + 6ujn 1
4ujn 11 + ujn 1
2
(86)
where  is an adjustable, non-negative parameter. Note that
ujn+21 4ujn+11 + 6ujn 1 4ujn 11 + ujn 21 = 4x 4(uxxxx)jn 1 + O(h6 )
= 4x 4(uxxxx)nj + O(h5 ) = O(h4 )
so that the term which is added, does not change the leading order truncation error, which in the form we
have written the equation, is O( 4t 3 ) = O(h3 ) (local/one-step truncation error).
A Von Neumann analysis of the modied scheme shows that, in addtion to the CFL condition   1, we
must have  < 1 for stability, and, further, that the per-step amplication factor for a mode with wave
number  is, to leading order
1  sin4 2
Thus the addition of the dissipation term is analagous to the use of an explicit \high frequency lter"
(low-pass lter), which has a fairly sharp rollover as  ! .
We note that an advantage to the use of explicit dissipation techniques (versus, for example, the use of an
intrinsically dissipative scheme) is that the amount of dissipation can be controlled by tuning the dissipation
parameter.

References
[1] Mitchell, A. R., and D. F. Griths, The Finite Dierence Method in Partial Dierential Equa-
tions, New York: Wiley (1980)
[2] Richtmeyer, R. D., and Morton, K. W., Dierence Methods for Initial-Value Problems, New
York: Interscience (1967)
[3] H.-O. Kreiss and J. Oliger, Methods for the Approximate Solution of Time Dependent Prob-
lems, GARP Publications Series No. 10, (1973)
[4] Gustatsson, B., H. Kreiss and J. Oliger, Time-Dependent Problems and Dierence Methods,
New York: Wiley (1995)
[5] Richardson, L. F., \The Approximate Arithmetical Solution by Finite Dierences of Physical Problems
involving Dierential Equations, with an Application to the Stresses in a Masonry Dam", Phil. Trans.
Roy. Soc., 210, (1910) 307{357.
[6] Anderson, E. em et al \Lapack Users' Guide - Release 2.0", (1994)
http://www.netlib.org/lapack/lug/lapack lug.html

23