02 Egr537-Lctrs
02 Egr537-Lctrs
02 Egr537-Lctrs
1
)
x
(
]
)
(
)
[J
(m+1
x =x
(m)
(m)
D f
f
0
i = i+1 f
i1
2h
LECTURES IN BASIC
COMPUTATIONAL
NUMERICAL ANALYSIS
J. M. McDonough
University of Kentucky
Lexington, KY 40506
E-mail: [email protected]
(
y =
y,t)
LECTURES IN BASIC
COMPUTATIONAL
NUMERICAL ANALYSIS
J. M. McDonough
Departments of Mechanical Engineering and Mathematics
University of Kentucky
c
1984,
1990, 1995, 2001, 2004, 2007
Contents
1 Numerical Linear Algebra
1.1 Some Basic Facts from Linear Algebra . . . . . . . . . . . . . . .
1.2 Solution of Linear Systems . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Numerical solution of linear systems: direct elimination .
1.2.2 Numerical solution of linear systems: iterative methods .
1.2.3 Summary of methods for solving linear systems . . . . . .
1.3 The Algebraic Eigenvalue Problem . . . . . . . . . . . . . . . . .
1.3.1 The power method . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Inverse iteration with Rayleigh quotient shifts . . . . . . .
1.3.3 The QR algorithm . . . . . . . . . . . . . . . . . . . . . .
1.3.4 Summary of methods for the algebraic eigenvalue problem
1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
1
5
5
16
24
25
26
29
30
31
32
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
33
33
33
34
38
38
39
41
41
43
43
3 Approximation Theory
3.1 Approximation of Functions . . . . . . . . . . . .
3.1.1 The method of least squares . . . . . . . .
3.1.2 Lagrange interpolation polynomials . . . .
3.1.3 Cubic spline interpolation . . . . . . . . .
3.1.4 Extraplotation . . . . . . . . . . . . . . .
3.2 Numerical Quadrature . . . . . . . . . . . . . . .
3.2.1 Basic NewtonCotes quadrature formulas
3.2.2 GaussLegendre quadrature . . . . . . . .
3.2.3 Evaluation of multiple integrals . . . . . .
3.3 Finite-Difference Approximations . . . . . . . . .
3.3.1 Basic concepts . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
45
45
45
49
55
60
60
60
65
66
68
68
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
ii
CONTENTS
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
68
70
71
72
73
76
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
77
77
77
80
90
94
99
101
102
103
104
108
109
110
114
118
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
119
119
120
121
122
124
124
128
130
136
144
145
148
149
149
155
159
3.4
3.5
3.6
References
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
160
List of Figures
1.1
1.2
1.3
1.4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
18
19
2.1
2.2
2.3
2.4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
36
37
39
40
3.1
3.2
3.3
3.4
3.5
3.6
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
46
53
54
55
61
65
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
= u
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
83
84
86
91
100
101
103
104
107
5.1
Methods for spatial discretization of partial differential equations; (a) finite difference, (b) finite element and (c) spectral. . . . . . . . . . . . . . . . . . . . . . . . .
Mesh star for forward-Euler method applied to heat equation . . . . . . . . . . . .
Matrix structure for 2-D CrankNicolson method. . . . . . . . . . . . . . . . . . .
Implementation of PeacemanRachford ADI. . . . . . . . . . . . . . . . . . . . . .
Matrix structure of centered discretization of Poisson/Dirichlet problem. . . . . . .
Analytical domain of dependence for the point (x, t). . . . . . . . . . . . . . . . . .
Numerical domain of dependence of the grid point (m, n + 1). . . . . . . . . . . . .
Difference approximation satisfying CFL condition. . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
122
125
138
141
146
150
152
156
5.2
5.3
5.4
5.5
5.6
5.7
5.8
iii
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Chapter 1
1.1
Before beginning our treatment of numerical solution of linear systems we will review a few important facts from linear algebra, itself. We typically think of linear algebra as being associated
with vectors and matrices in some finite-dimensional space. But, in fact, most of our ideas extend
quite naturally to the infinite-dimensional spaces frequently encountered in the study of partial
differential equations.
We begin with the basic notion of linearity which is crucial to much of mathematical analysis.
Definition 1.1 Let S be a vector space defined on the real numbers R (or the complex numbers
C), and let L be an operator (or transformation) whose domain is S. Suppose for any u, v S and
a, b R (or C) we have
L(au + bv) = aLu + bLv.
(1.1)
Then L is said to be a linear operator.
Examples of linear operators include M N matrices, differential operators and integral operators.
It is generally important to be able to distinguish linear and nonlinear operators because problems involving only the former can often be solved without recourse to iterative procedures. This
1
is seldom true for nonlinear problems, with the consequence that corresponding algorithms must
be more elaborate. This will become apparent as we proceed.
One of the most fundamental properties of any object, be it mathematical or physical, is its
size. Of course, in numerical analysis we are always concerned with the size of the error in any
particular numerical approximation, or computational procedure. There is a general mathematical
object, called the norm, by which we can assign a number corresponding to the size of various
mathematical entities.
Definition 1.2 Let S be a (finite- or infinite-dimensional) vector space, and let k k denote the
mapping S R+ {0} with the following properties:
i) kvk 0, v S with kvk = 0 iff v 0,
ii) kavk = |a| kvk,
v S, a R,
v, w S.
kvk2 =
N
X
vi2
i=1
!12
(1.2)
The proof that k k2 , often called the Euclidean norm, or simply the 2-norm, satisfies the three
conditions of the definition is straightforward, and is left to the reader. (We note here that it is
common in numerical analysis to employ the subscript E to denote this norm and use the subscript
2 for the spectral norm of matrices. But we have chosen to defer to notation more consistent
with pure mathematics.) Another useful norm that we often encounter in practice is the max norm
or infinity norm defined as
kvk = max |vi | .
(1.3)
1iN
In the case of Euclidean spaces, we can define another useful object related to the Euclidean
norm, the inner product (often called the dot product when applied to finite-dimensional vectors).
Definition 1.3 Let S be a N-dimensional Euclidean space with v, w S. Then
hv, wi
N
X
vi wi
(1.4)
i=1
Theorem 1.1 (CauchySchwarz) Let S be an inner-product space with inner product h , i and
norm k k2 . If v, w S, then
hv, wi kvk2 kwk2 .
(1.5)
We have thus far introduced the 2-norm, the infinity norm and the inner product for spaces of
finite-dimensional vectors. It is worth mentioning that similar definitions hold as well for infinitedimensional spaces, i.e., spaces of functions. For example, suppose f (x) is a function continuous
on the closed interval [a, b], denoted f C[a, b]. Then
kf k = max |f (x)|.
(1.6)
x [a,b]
Z
f dx
12
The space consisting of all functions f such that kf k2 < is the canonical Hilbert space, L2 [a, b].
The CauchySchwarz inequality holds in any such space, and takes the form
Z
b
a
f g dx
Z
f dx
a
12 Z
g dx
a
21
f, g L2 [a, b].
We next need to consider some corresponding ideas regarding specific calculations for norms of
operators. The general definition of an operator norm is as follows.
Definition 1.4 Let A be an operator whose domain is D. Then the norm of A is defined as
kAk max kAxk .
(1.7)
kxk=1
xD(A)
kxk6=0
xD(A)
kAxk
,
kxk
from which follows an inequality similar to the CauchySchwarz inequality for vectors,
kAxk kAkkxk.
(1.8)
We should remark here that (1.8) actually holds only in the case when the matrix and vector norms
appearing in the expression are compatible, and this relationship is often used as the definition
of compatibility. We will seldom need to employ this concept in the present lectures, and the reader
is referred to, e.g., Isaacson and Keller [15] (Chap. 1) for additional information.
We observe that neither (1.7) nor the expression following it is suitable for practical calculations;
we now present three norms that are readily computed, at least for M N matrices. The first of
these is the 2-norm, given in the matrix case by
kAk2
M,N
X
i,j=1
a2ij
(1.9)
Two other norms are also frequently employed. These are the 1-norm
kAk1 = max
1jN
1iM
M
X
i=1
|aij | ,
N
X
j=1
|aij | .
(1.10)
(1.11)
We note that although the definition of the operator norm given above was not necessarily finitedimensional, we have here given only finite-dimensional practical computational formulas. We will
see later that this is not really a serious restriction because problems involving differential operators,
one of the main instances where norms of infinite-dimensional operators are needed, are essentially
always solved via discrete approximations leading to finite-dimensional matrix representations.
There is a final, general comment that should be made regarding norms. It arises from the
fact, mentioned earlier, that in any given vector space many different norms might be employed.
A comparison of the formulas in Eqs. (1.2) and (1.3), for example, will show that the number one
obtains to quantify the size of a mathematical object, a vector in this case, will change according to
which formula is used. Thus, a reasonable question is, How do we decide which norm to use? It
turns out, for the finite-dimensional spaces we will deal with herein, that it really does not matter
which norm is used, provided only that the same one is used when making comparisons between
similar mathematical objects. This is the content of what is known as the norm equivalence theorem:
all norms are equivalent on finite-dimensional spaces in the sense that if a sequence converges in
one norm, it will converge in any other norm (see Ref. [15], Chap. 1). This implies that in practice
we should usually employ the norm that requires the least amount of floating-point arithmetic
for its evaluation. But we note here that the situation is rather different for infinite-dimensional
spaces. In particular, for problems involving differential equations, determination of the function
space in which a solution exists (and hence, the appropriate norm) is a significant part of the overall
problem.
We will close this subsection on basic linear algebra with a statement of the problem whose
numerical solution will concern us throughout most of the remainder of this chapter, and provide
the formal, exact solution. We will study solution procedures for the linear system
Ax = b,
(1.12)
(1.13)
It was apparently not clear in the early days of numerical computation that direct application
of (1.13), i.e., computing A1 and multiplying b, was very inefficientand this approach is rather
natural. But if A is a NN matrix, as much as O(N 4 ) floating-point arithmetic operations may be
required to produce A1 . On the other hand, if the Gaussian elimination procedure to be described
in the next section is used, the system (1.12) can be solved for x, directly, in O(N 3 ) arithmetic
operations. In fact, a more cleverly constructed matrix inversion routine would use this approach
to obtain A1 in O(N 3 ) arithmetic operations, although the precise number would be considerably
greater than that required to directly solve the system. It should be clear from this that one should
never invert a matrix to solve a linear system unless the inverse matrix, itself, is needed for other
purposes, which is not usually the case for the types of problems treated in these lectures.
1.2
In this section we treat the two main classes of methods for solving linear systems: i) direct
elimination, and ii) iterative techniques. For the first of these, we will consider the general case
of a nonsparse N N system matrix, and then study a very efficient elimination method designed
specifically for the solution of systems whose matrices are sparse, and banded. The study of the
second topic, iterative methods, will include only very classical material. It is the authors opinion
that students must be familiar with this before going on to study the more modern, and much
more efficient, methods. Thus, our attention here will be restricted to the topics Jacobi iteration,
GaussSeidel iteration and successive overrelaxation.
1.2.1
a11
a21
a31
(1.14)
where the matrix A with elements aij is assumed to be nonsingular. Moreover, we assume that
no aij or bi is zero. (This is simply to maintain complete generality.) If we perform the indicated
matrix/vector multiplication on the left-hand side of (1.14) we obtain
a11 x1 + a12 x2 + a13 x3 = b1 ,
a21 x1 + a22 x2 + a23 x3 = b2 ,
(1.15)
b3
,
a33
and working backward, in order, to x1 . This motivates trying to find combinations of the equations
in (1.15) such that the lower triangle of the matrix A is reduced to zero. We will see that the
resulting formal procedure, known as Gaussian elimination, or simply direct elimination, is nothing
more than a systematic approach to methods from high school algebra, organized to lend itself to
machine computation.
The Gaussian elimination algorithm proceeds in stages. At each stage a pivot element is selected
so as to appear in the upper left-hand corner of the largest submatrix not processed in an earlier
stage. Then all elements in the same column below this pivot element are zeroed by means of
elementary operations (multiplying by a constant, replacing an equation with a linear combination
of equations). This continues until the largest submatrix consists of only the element in the lower
right-hand corner. We demonstrate this process with the system (1.15).
At stage 1, the pivot element is a11 , and in carrying out calculations for this stage we must zero
the elements a21 and a31 . If we multiply the first equation by the Gauss multiplier
m21
a21
,
a11
a21 a11
a21
a11
x1 + a22 a12
a21
a11
a21
a21
x2 + a23 a13
x3 = b 2 b 1
,
a11
a11
or
0 x1 + (a22 m21 a12 )x2 + (a23 m21 a13 )x3 = b2 m21 b1 .
We see that if we replace the second equation of (1.15) with this result we will have achieved the
goal of zeroing a21 .
Similarly, define m31 a31 /a11 , and carry out analogous operations between the first and third
equations of (1.15). This completes the first stage, because we have now zeroed all elements below
the pivot element in the first column of the matrix A. The system can now be expressed as
a11 a12 a13
x1
b1
0 a22 a23 x2 = b2 ,
0 a32 a33
x3
b3
where
b2 = b2 m21 b1 ,
b3 = b3 m31 b1 .
Itis worthwhile to observe that no element of the first row (equation) has been changed by this
process; but, in general, all other elements will have been changed.
We are now prepared to continue to stage 2. Now the pivot element is a22 (not a22 !) and we
have only to zero the element a32 . To show more clearly what is happening, we write the part of
the system being considered at this stage as
a22 x2 + a23 x3 = b2 ,
a32 x2 + a33 x3 = b3 .
Now we define m32 a32 /a22 , multiply this by the first equation (of the current stage), and subtract
from the second. This yields
a32
a32
a32
x2 + a33 a23
x3 = b3 b2
.
a32 a22
a22
a22
a22
Thus, if we define
a
33 a33 m32 a23 ,
and
b
3 b3 m32 b2 ,
a11 a12 a13
x1
b1
0 a22 a23 x2 = b2 .
0
0 a
x3
b
33
3
(1.16)
In summarizing this elimination process there are several important observations to be made
regarding the form of (1.16), and how it was achieved. First, we note that the matrix in (1.16) is
of the form called upper triangular; and as suggested earlier, such a system can be easily solved. In
arriving at this form, the first equation in the original system was never changed, as already noted,
but more importantly, neither was the solution vector, or the ordering of its elements, changed.
Thus, the solution to (1.16) is identical to that of (1.14). Obviously, if this were not true, this whole
development would be essentially worthless. Next we observe that for the 33 matrix treated here,
only two elimination stages were needed. It is easily seen that this generalizes to arbitrary N N
matrices for which N 1 elimination stages are required to achieve upper triangular form. Finally,
we point out the tacit assumption that the pivot element at each stage is nonzero, for otherwise the
Gauss multipliers, mij would not be defined. It is quite possible for this assumption to be violated,
even when the matrix A is nonsingular. We will treat this case later.
We can now complete the solution of Eq. (1.14) via the upper triangular system (1.16). We
have
x3 = b
3 /a33 ,
Do i = k + 1, N
mik = aik /akk
bi = bi mik bk
Do j = k + 1, N
aij = aij mik akj
Repeat j
xi = 0
Do j = i + 1, N
xi = xi + aij xj
Repeat j
xi = (bi xi )/aii
Repeat i
At this point we comment that it should be clear from the structure of this algorithm that
O(N 3 ) arithmetic operations are required to obtain a solution to a linear system using Gaussian
elimination. In particular, in the forward elimination step there is a nesting of three DO loops, each
of which runs O(N ) times. In addition, the backward substitution step requires O(N 2 ) operations;
but for large N this is negligible compared with O(N 3 ). It is important to realize that even on
modern supercomputers, there are many situations in which this amount of arithmetic is prohibitive,
and we will investigate ways to improve on this operation count. However, for nonsparse matrices
Gaussian elimination is in most cases the preferred form of solution procedure.
We now state a theorem that provides conditions under which Gaussian elimination is guaranteed to succeed.
Theorem 1.2 (LU Decomposition) Given a square NN matrix A, let Ak denote the principal minor constructed from the first k rows and k columns. Suppose that det(Ak ) 6= 0 k = 1, 2, . . . , N 1.
Then a unique lower triangular matrix L = (ij ) with 11 = 22 = = N N = 1 and a unique
upper triangular matrix U = (uij ) LU = A. Furthermore,
det(A) =
N
Y
uii .
i=1
We remark that it is not entirely obvious that Gaussian elimination, as presented above, actually
corresponds to LU decomposition of the theorem; but, in fact, the lower triangular matrix of the
theorem can be recovered from the Gauss multipliers, as demonstrated in Forsythe and Moler [9],
for example. The importance of LU decomposition, per se, will be more apparent later when we
treat direct solution of sparse systems having band matrices.
Rounding Errors in Gaussian Elimination
We next consider a well-known example (see Ref. [9]), the purpose of which is to demonstrate the
effects of rounding errors in Gaussian elimination. The system
1.0
0.0001 1.0 x1
=
(1.17)
2.0
1.0
1.0 x2
has the solution x1 = 1.000, x2 = 1.000 when solved by Cramers rule with rounding to three
significant digits. We now perform Gaussian elimination with the same arithmetic precision. Since
this is a 22 system, there will be only a single forward elimination stage. We choose the natural
pivot element, a11 = 0.0001, and compute the Gauss multiplier
m21 = a21 /a11 = 104 .
Then, to three significant digits,
a22 = a22 m21 a12 = 1.0 (104 )(1.0) = 104 ,
and
b2 = b2 m21 b1 = 2.0 (104 )(1.0) = 104 .
We then have
x2 = b2 /a22 = 1.0,
and
1.0 (1.0)(1.0)
b1 a12 x2
=
= 0.0.
a11
104
Thus, the second component of the solution is accurate, but the first component is completely
wrong.
An obvious remedy is to increase the precision of the arithmetic, since three significant digits
is not very accurate to begin with. In particular, the rounding errors arise in this case during
subtraction of a large number from a small one, and information from the small number is lost
due to low arithmetic precision. Often, use of higher precision might be the simplest approach
to handle this. Unfortunately, even for the simple 2 2 system just treated we could choose a11
in such a way that Gaussian elimination would fail on any finite machine. The desire to produce
machine-independent algorithms led numerical analysts to study the causes of rounding errors
in direct elimination methods. It is not too surprising that the error originates in the forward
elimination step, and is caused by relatively (with respect to machine precision) small values of the
pivot element used at any particular stage of the elimination procedure which, in turn, can result
in large Gauss multipliers, as seen in the preceding example. From this piece of information it is
clear that at the beginning of each stage it would be desirable to arrange the equations so that the
pivot element for that stage is relatively large.
There are three main strategies for incorporating this idea: i) column pivoting, ii) row pivoting
and iii) complete pivoting. The first of these is the easiest to implement because it requires only row
interchanges, and therefore does not result in a re-ordering of the solution vector. Row pivoting
requires column interchanges, and thus does reorder the solution vector; as would be expected,
complete pivoting utilizes both row and column interchanges.
Numerical experiments, performed for a large collection of problems, have shown that complete
pivoting is more effective than either of the partial pivoting strategies, which are generally of
comparable effectiveness. Thus, we recommend row interchange (column pivoting) if a partial
pivoting strategy, which usually is sufficient, is employed. The Gaussian elimination algorithm
given earlier requires the following additional steps, inserted between the k and i DO loops of
the forward elimination step, in order to implement the row interchange strategy.
x1 =
10
It is important to note that these steps must be executed at every stage of the forward elimination procedure. It is not sufficient to seek the largest pivot element prior to the first stage, and
then simply use the natural pivot elements thereafter. It should be clear from the detailed analysis
of the 3 3 system carried out earlier that the natural pivot element at any given stage may be
very different from the element in that position of the original matrix.
We should also point out that there are numerous other procedures for controlling round-off
errors in Gaussian elimination, treatment of which is beyond the intended scope of these lectures.
Some of these, such as balancing of the system matrix, are presented in the text by Johnson and
Riess [16], for example.
Condition Number of a Matrix
We shall conclude this treatment of basic Gaussian elimination with the derivation of a quantity
which, when it can be calculated, provides an a priori estimate of the effects of rounding and data
errors on the accuracy of the solution to a given linear system. This quantity is called the condition
number of the system matrix.
We begin with the general system
Ax = b
with exact solution x = A1 b. We suppose that A can be represented exactly, and that somehow
A1 is known exactly. We then inquire into the error in the solution vector, x, if the forcing vector
b is in error by b. If the error in x due to b is denoted x, we have
A(x + x) = b + b,
from which it follows that
x = A1 b.
Thus, (assuming use of compatible matrix and vector normsrecall Eq. (1.8))
kxk kA1 k kbk .
(1.18)
11
(1.19)
The inequality in (1.19) shows that the relative error in the data vector b is amplified by an
amount kAkkA1 k to produce a corresponding relative error in the solution vector. This quantity
is called the condition number of A, denoted
cond(A) kAk kA1 k ,
(1.20)
or often simply (A). With a bit more work, it can also be shown (see [9], for example) that
kAk
kxk
cond(A)
.
kx + xk
kAk
(1.21)
This inequality, together with (1.19), shows the gross effects of rounding errors on the numerical
solution of linear systems. In this regard it is of interest to note that the elimination step of Gaussian
elimination results in the replacement of the original system matrix A with the LU decomposition
LU A, and the original forcing vector b with b b, where A and b are due to round-off errors.
There are several remarks to make regarding the condition number. The first is that its value
depends on the norm in which it is calculated. Second, determination of its exact value depends on
having an accurate A1 . But if A is badly conditioned, it will be difficult (and costly) to compute
A1 accurately. Hence, cond(A) will be inaccurate. There has been considerable research into
obtaining accurate and efficiently-computed approximations to cond(A). (The reader is referred to
Dongarra et al. [6] for computational procedures.) A readily calculated rough approximation is the
following:
kAk2
,
(1.22)
cond(A)
| det(A)|
which is given by Hornbeck [14] in a slightly different context. It should be noted, however, that
this is only an approximation, and not actually a condition number, since it does not satisfy the
easily proven inequality,
cond(A) 1 .
A system is considered badly conditioned, or ill conditioned, when cond(A) 1, and well conditioned when cond(A) O(1). However, ill conditioning must be viewed relative to the arithmetic
precision of the computer hardware being employed, and relative to the required precision of computed results. For example, in current single-precision arithmetic (32-bit words), condition numbers
up to O(103 ) can typically be tolerated; in double precision (64-bit words) accurate solutions can
be obtained even when cond(A) & O(106 ). The interested reader may find further discussions of
this topic by Johnson and Riess [16] quite useful.
12
of (mostly) nonzero elements, and all other elements are zero. This particular matrix structure
is still too general to permit application of direct band elimination. In general, this requires that
matrices be in the form called compactly banded, shown in Fig. 1.2. Part (a) of this figure shows
(a)
(b)
a three-band, or tridiagonal matrix, while part (b) displays a pentadiagonal matrix. Both of these
frequently occur in practice in the contexts mentioned above; but the former is more prevalent, so
it will be the only one specifically treated here. It will be clear, however, that our approach directly
13
extends to more general cases. In fact, it is possible to construct a single algorithm that can handle
any compactly-banded matrix.
The approach we shall use for solving linear systems with coefficient matrices of the above form
is formal LU decomposition as described in Theorem 1.2. The algorithm we obtain is completely
equivalent (in terms of required arithmetic and numerical stability with respect to rounding errors)
to Gaussian elimination, and to the well-known Thomas algorithm [34] for tridiagonal systems. We
write the system as
a1,2 a1,3
x1
b1
x2
b2
x3
b3
.
.
.
.
xi
.
.
.
aN1,1 aN 1,2 aN1,3
aN,1
aN,2
bi
.
.
.
xN1
bN1
xN
bN
(1.23)
Observe that indexing of the elements of A is done differently than in the nonsparse case treated
earlier. In particular, the first index in the present sparse case is the equation (row) index, as was
true for nonsparse matrices; but the second index is now the band number, counting from lower
left to upper right corners of the matrix, rather than the column number.
In order to derive the LU decomposition formulas, we assume the matrix A can be decomposed
into the lower- and upper-triangular matrices, L and U , respectively (as guaranteed by the LUdecomposition Theorem). We then formally multiply these matrices back together, and match
elements of the resulting product matrix LU with those of the original matrix A given in Eq.
(1.23). This will lead to expressions for the elements of L and U in terms of those of A.
A key assumption in this development is that L and U are also band matrices with the structure
of the lower and upper, respectively, triangles of the original matrix A. We thus have
1
d1
e1
1
c2 d 2
e2
1
c3 d3
L=
e3
, U=
dN1
cN dN
.
1
eN1
1
(1.24)
14
d1
d1 e1
c2 d2+ c 2 e1 d2e2
c3
LU =
d3+ c3 e2 d3 e3
ci
di + ci ei1 di ei
dN1 eN1
cN
dN+ cN eN1
Thus, as we would hope, the matrices of our proposed LU decomposition when multiplied
together lead to the same tridiagonal matrix structure as that of the original matrix A in Eq.
(1.23), and from the LU-decomposition theorem we know that the matrix elements must equal
those of the original matix Aelement by element. Thus, we immediately see that ci = ai,1 must
hold i = 2, . . . , N . Also, d1 = a12 , which permits calculation of e1 = a13 /d1 = a13 /a12 . In
general, we see that we can immediately calculate di and ei at each row of LU , since ci and ei1
are already known. In particular, we have
di = ai,2 ci ei1
and
ei1 = ai1,3 /di1 .
Hence, determination of the components of L and U is a straightforward, well-defined process.
Once we have obtained an LU decomposition of any matrix, whether or not it is sparse or
banded, solution of the corresponding system is especially simple. Consider again the general
system
Ax = b,
and suppose we have obtained upper and lower triangular matrices as in the LU-decomposition
Theorem, i.e., such that A = LU . Then
LU x = b,
and if we set y = U x, we have
Ly = b.
But b is given data, and L is a lower triangular matrix; hence, we can directly solve for y by forward
substitution, starting with the first component of y. Once y has been determined, we can find the
desired solution, x, from U x = y via backward substitution, since U is upper triangular. Again,
we emphasize that this procedure works for any LU decomposition; its use is not restricted to the
tridiagonal case that we have treated here.
We now summarize the preceding derivation for solution of tridiagonal systems with the following pseudo-language algorithm.
15
xi = bi ai,3 xi+1
Repeat i
It is worthwhile to compare some of the details of this algorithm with Gaussian elimination
discussed earlier for nonsparse systems. Recall in that case we found, for a N N system matrix,
that the total required arithmetic was O(N 3 ), and the necessary storage was O(N 2 ). For tridiagonal
systems the situation is much more favorable. Again, for a N N system (i.e., N equations in N
unknowns) only O(N ) arithmetic operations are needed to complete a solution. The reason for
this significant reduction of arithmetic can easily be seen by noting that in the algorithm for the
tridiagonal case there are no nested DO loops. Furthermore, the longest loop is only of length
N . It is clear from (1.23) that at most 5N words of storage are needed, but in fact, this can be
reduced to 4N 2 by storing the solution vector, x, in the same locations that originally held the
right-hand side vector, b. Rather typical values of N are in the range O(102 ) to O(103 ). Hence, very
significant savings in arithmetic (and storage) can be realized by using band elimination in place of
the usual Gaussian elimination algorithm. We should note, however, that it is rather difficult (and
seldom done) to implement pivoting strategies that preserve sparsity. Consequently, tridiagonal LU
decomposition is typically applied only in those situations where pivoting would not be necessary.
A final remark, of particular importance with respect to modern vector and parallel supercomputer architectures, should be made regarding the tridiagonal LU-decomposition algorithm just
presented. It is that this algorithm can be neither vectorized nor parallelized in any direct manner;
it works best for scalar processing. There are other band-matrix solution techniques, generally
known as cyclic reduction methods, that can be both vectorized and parallelized to some extent.
Discussion of these is beyond the intended scope of the present lectures; the reader is referred to
Duff et al. [8], for example, for more information on this important topic. On the other hand,
parallelization of application of the algorithm (rather than of the algorithm, itself) is very effective,
and widely used, in numerically solving partial differential equations by some methods to be treated
in Chap. 5 of these notes.
16
1.2.2
From the preceding discussions of direct elimination procedures it is easily seen that if for any
reason a system must be treated as nonsparse, the arithmetic required for a solution can become
prohibitive already for relatively small N , say N O(105 ). In fact, there are many types of
sparse systems whose structure is such that the sparsity cannot be exploited in a direct elimination
algorithm. Figure 1.1 depicts an example of such a matrix. Matrices of this form arise in some finitedifference and finite-element discretizations of fairly general two-dimensional elliptic operators, as
we shall see in Chap. 5; for accurate solutions one may wish to use at least 104 points in the
finite-difference mesh, and it is now not unusual to employ as many as 106 points, or even more.
This results in often unacceptably long execution times, even on supercomputers, when employing
direct elimination methods. Thus, we are forced to attempt the solution of this type of problem
via iterative techniques. This will be the topic of the present section.
We will begin with a general discussion of fixed-point iteration, the basis of many commonly
used iteration schemes. We will then apply fixed-point theory to linear systems, first in the form
of Jacobi iteration, then with an improved version of this known as GaussSeidel iteration, and
finally with the widely-used successive overrelaxation.
Fundamentals of Fixed-Point Iteration
Rather than start immediately with iterative procedures designed specifically for linear problems, as
is usually done in the present context, we will begin with the general underlying principles because
these will be of use later when we study nonlinear systems, and because they provide a more direct
approach to iteration schemes for linear systems as well. Therefore, we will briefly digress to study
some elementary parts of the theory of fixed points.
We must first introduce some mathematical notions. We generally view iteration schemes as
methods which somehow generate a sequence of approximations to the desired solution for a given
problem. This is often referred to as successive approximation, or sometimes as trial and error.
In the hands of a novice, the latter description may, unfortunately, be accurate. Brute-force, trialand-error methods often used by engineers are almost never very efficient, and very often do not
work at all. They are usually constructed on the basis of intuitive perceptions regarding the physical
system or phenomenon under study, with little or no concern for the mathematical structure of the
equations being solved. This invites disaster.
In successive approximation methods we start with a function, called the iteration function,
which maps one approximation into another, hopefully better, one. In this way a sequence of possible solutions to the problem is generated. The obvious practical question that arises is, When have
we produced a sufficient number of approximations to have obtained an acceptably accurate answer? This is simply a restatement of the basic convergence question from mathematical analysis,
so we present a few definitions and notations regarding this.
N
Definition 1.5 Let {ym }
m=1 be a sequence in R . The sequence is said to converge to the limit
y RN if > 0 M (depending on ) m M , kyym k < . We denote this by lim ym = y .
m
We note here that the norm has not been specified, and we recall the earlier remark concerning
equivalence of norms in finite-dimensional spaces. (More information on this can be found, for
example, in Apostol [1].) We also observe that when N = 1, we merely replace norm k k with
absolute value | | .
It is fairly obvious that the above definition is not generally of practical value, for if we knew
the limit y, which is required to check convergence, we probably would not need to generate the
17
sequence in the first place. To circumvent such difficulties mathematicians invented the notion of
a Cauchy sequence given in the following definition.
N
Definition 1.6 Let {ym }
m=1 R , and suppose that > 0 M (depending on ) m, n M ,
kym yn k < . Then {ym } is a Cauchy sequence.
By itself, this definition would not be of much importance; but it is a fairly easily proven fact
from elementary analysis (see, e.g., [1]) that every Cauchy sequence in a complete metric space
converges to an element of that space. It is also easy to show that RN is a complete metric space
N < . Thus, we need only demonstrate that our successive approximations form a Cauchy
sequence, and we can then conclude that the sequence converges in the sense of the earlier definition.
Although this represents a considerable improvement over trying to use the basic definition
in convergence tests, it still leaves much to be desired. In particular, the definition of a Cauchy
sequence requires that kym yn k < hold > 0, and m, n M , where M , itself, is not
specified, a priori. For computational purposes it is completely unreasonable to choose smaller
than the absolute normalized precision (the machine ) of the floating-point arithmetic employed.
It is usually sufficient to use values of O(e/10), where e is the acceptable error for the computed
results. The more difficult part of the definition to satisfy is m, n M . However, for wellbehaved sequences, it is typically sufficient to choose n = m + k where k is a specified integer
between one and, say 100. (Often k = 1 is used.) The great majority of computer-implemented
iteration schemes test convergence in this way; that is, the computed sequence {ym } is considered
to be converged when kym+1 ym k < for some prescribed (and often completely fixed) . In many
practical calculations 103 represents quite sufficient accuracy, but of course this is problem
dependent.
Now that we have a means by which the convergence of sequences can be tested, we will study a
systematic method for generating these sequences in the context of solving equations. This method
is based on a very powerful and basic notion from mathematical analysis, the fixed point of a
function, or mapping.
(1.25)
and in this form we recognize that a fixed point of f is a zero (or root) of the function g(x) xf (x).
Hence, if we can find a way to compute fixed points, we automatically obtain a method for solving
equations.
From our earlier description of successive approximation, intuition suggests that we might try
18
xm = f (xm1 )
where x0 is an initial guess. This procedure generates the sequence {xm } of approximations to the
fixed point x , and we continue this until kxm+1 xm k < .
It is of interest to analyze this scheme graphically for f : D D, D R1 . This is presented in
the sketch shown as Fig. 1.3. We view the left- and right-hand sides of x = f (x) as two separate
functions, y = x and z = f (x). Clearly, there exists x = x such that y = z is the fixed point of
f ; that is, x = f (x ). Starting at the initial guess x0 , we find f (x0 ) on the curve z = f (x). But
according to the iteration scheme, x1 = f (x0 ), so we move horizontally to the curve y = x. This
locates the next iterate, x1 on the x-axis, as shown in the figure. We now repeat the process by
again moving vertically to z = f (x) to evaluate f (x1 ), and continuing as before. It is easy to see
that the iterations are very rapidly converging to x for this case.
f ( x)
y=x
z = f ( x)
f (x* )
f ( x 1)
x* = f (x*)
f ( x 0)
x0
x1
x 2 x*
19
in Fig. 1.4. We first observe that the function f indeed has a fixed point, x . But the iterations
starting from x0 < x do not converge to this point. Moreover, the reader may check that the
iterations diverge also for any initial guess x0 > x . Comparison of Figs. 1.3 and 1.4 shows that
in Fig. 1.3, f has a slope less than unity (in magnitude) throughout a neighborhood containing
the fixed point, while this is not true for the function in Fig. 1.4. An iteration function having a
z = f (x )
f ( x)
y=x
f ( x*)
x* = f (x* )
f (x0)
f (x1)
x2
x1
x0 x*
x, y D .
(1.26)
x, y D ,
is called a Lipschitz condition, and L is the Lipschitz constant. Any function f satisfying such a
condition is said to be a Lipschitz function.
20
There are several things to note regarding the above theorem. The first is that satisfaction of
the Lipschitz condition with L < 1 is sufficient, but not always necessary, for convergence of the
corresponding iterations. Second, for any set D, and mapping f with (1.26) holding throughout,
x is the unique fixed point of f in D. Furthermore, the iterations will converge to x using any
starting guess, whatever, so long as it is an element of the set D. Finally, the hypothesis that f is
continuous in D is essential. It is easy to construct examples of iteration functions satisfying all of
the stated conditions except continuity, and for which the iterations fail to converge.
Clearly, it would be useful to be able to calculate the Lipschitz constant L for any given function
f . This is not always possible, in general; however, there are some practical situations in which L
can be calculated. In particular, the theorem requires only that f be continuous on D, but if we
assume further that f possesses a bounded derivative in this domain, then the mean value theorem
gives the following for D R1 :
f (b) f (a) = f ()(b a) ,
Then
max |f (x)| |b a|,
x[a,b]
and we take
L = max |f (x)|.
x[a,b]
(We comment that, for technical reasons associated with the definition of the derivative, we should
take D = [a , b + ], > 0, in this case.)
For D RN , the basic idea is still the same, but more involved. The ordinary derivative must be
replaced by a Frechet derivative, which in the finite-dimensional case is just the Jacobian matrix,
and the max operation must be applied to the largest (in absolute value) eigenvalue of this matrix,
viewed as a function of x D RN . We note that there are other possibilities for computing L in
multi-dimensional spaces, but we will not pursue these here.
Jacobi Iteration
We are now prepared to consider the iterative solution of systems of linear equations. Thus, we
again study the equation Ax = b. For purposes of demonstration we will consider the same general
33 system treated earlier in the context of Gaussian elimination:
a11 x1 + a12 x2 + a13 x3 = b1
a21 x1 + a22 x2 + a23 x3 = b2
(1.27)
21
It is clear that with x = (x1 , x2 , x3 )T and f = (f1 , f2 , f3 )T , we have obtained the form x = f (x),
as desired. Moreover, the iteration scheme corresponding to successive approximation can be
expressed as
1
(m+1)
(m)
(m)
x1
=
b1 a12 x2 a13 x3
a11
1
(m)
(m)
(m+1)
b2 a21 x1 a23 x3
x2
=
(1.28)
a22
1
(m)
(m)
(m+1)
b3 a31 x1 a32 x2
x3
=
,
a33
(m+1)
xi
1
aii
N
X
(m)
bi
aij xj
fi (x1 , x2 , . . . , xN ).
(1.29)
j=1
j6=i
xi
Repeat i
1
aii
N
X
(m)
bi
aij xj
j=1
j6=i
3. Test convergence:
(m+1)
(m)
xi < , then stop
If max xi
i
(m)
else set xi
(m+1)
= xi
, i = 1, . . . , N
m=m+1
go to 2
N
X
j=1
j6=i
|aij |
i = 1, 2, . . . , N.
(1.30)
22
The condition given by inequality (1.30) is known as (strict) diagonal dominance of the matrix
A. It has not been our previous practice to prove theorems that have been quoted. However, the
proof of the above theorem provides a clear illustration of application of the contraction mapping
principle, so we will present it here.
Proof. We first observe that any square matrix A can be decomposed as A = D L U , where
D is a diagonal matrix with elements equal to those of the main diagonal of A, and L and U are
respectively lower and upper triangular matrices with zero diagonals, and whose (nonzero) elements
are the negatives of the corresponding elements of A. (Note that these matrices are not the L and
U of the LU-decomposition theorem stated earlier.) Then (D L U )x = b, and
(1.31)
x(m+1) = D 1 (L + U )x(m) + D 1 b f x(m) .
It is easily checked that this fixed-point representation is exactly the form of (1.28) and (1.29).
To guarantee convergence of the iteration scheme (1.31) we must find conditions under which
the corresponding Lipschitz constant has magnitude less than unity. From (1.31) it follows that
the Lipschitz condition is
kf (x) f (y)k =
D 1 (L + U )x + D 1 b D 1 (L + U )y D 1 b
=
D 1 (L + U )(x y)
D 1 (L + U )
kx yk .
K =
D 1 (L + U )
,
N
1 X
K = max
|aij |
.
i
|aii |
j=1
j6=i
N
X
j=1
j6=i
|aij | i = 1, 2, . . . , N.
23
through all equations of the system. It turns out that this is, in fact, usually true, and it provides
the basis for the method known as GaussSeidel iteration. In this scheme the general ith equation
is given by
N
i1
X
X
1
(m)
(m+1)
(m+1)
aij xj ,
i = 1, 2, . . . , N .
(1.32)
aij xj
bi
xi
=
aii
j=i+1
j=1
We note that this is not strictly a fixed-point algorithm because both x(m) and x(m+1) appear in
the iteration function. It is fairly easy, however, to effect a rearrangement that formally achieves
the fixed-point form. This is important for theoretical studies but is of little value for practical
computation, so we will not pursue it here.
It is worthwhile to make some comparisons with Jacobi iteration. Because all equations are
updated simultaneously in Jacobi iteration, the rate of convergence (or divergence) is not influenced
by the order in which individual solution components are evaluated. This is not the case for Gauss
Seidel iterations. In particular, it is possible for one ordering of the equations to be convergent,
while a different ordering is divergent. A second point concerns the relationship between rates of
convergence for Jacobi and GaussSeidel iterations; it can be shown theoretically that when both
methods converge, GaussSeidel does so twice as fast as does Jacobi for typical problems having
matrix structure such as depicted in Fig. 1.1. This is observed to a high degree of consistency in
actual calculations, the implication being that one should usually employ GaussSeidel instead of
Jacobi iterations.
Successive Overrelaxation
One of the more widely-used methods for iterative solution of sparse systems of linear equations is
successive overrelaxation (SOR). This method is merely an accelerated version of the GaussSeidel
procedure discussed above. Suppose we let xi represent the ith component of the solution obtained
using GaussSeidel iteration, Eq. (1.32). Let be the so-called relaxation parameter. Then we can
often improve the convergence rate of the iterations by defining the weighted average
(m+1)
xi
(m)
= (1 )xi + xi
(m)
(m)
= xi + xi xi
(m)
= xi
+ xi .
Then, replacing xi with (1.32) leads to the general ith equation for SOR:
N
i1
X
X
1
(m)
(m)
(m+1)
(m+1)
(m)
aij xj xi ,
aij xj
xi
= x i + bi
aii
j=1
(1.33)
j=i+1
i =1, 2, . . . , N .
For 1 < < 2, this corresponds to SOR; for 0 < < 1, we obtain a damped form of Gauss
Seidel often termed underrelaxation. When writing a GaussSeidel code one should always include
the relaxation parameter as an input. Clearly Eq. (1.33), SOR, reduces to (1.32), GaussSeidel,
when = 1. Thus, a computer program written using (1.33) exhibits a great deal of flexibility in
terms of controlling the rate of convergence of the iterations.
We now present a pseudo-language algorithm for implementing SOR.
24
1
bi
aii
i1
X
(m+1)
aij xj
j=1
N
X
xi
(m)
= xi
(m)
aij xj
j=i+1
(m)
xi
+ xi
Repeat i
4. Test convergence
if maxdif < , then print results, and stop
Repeat m
5. Print message that iterations failed to converge in maxitr iterations
A few comments should be made regarding the above pseudo-language algorithm. First, it
should be noticed that we have provided a much more detailed treatment of the SOR algorithm in
comparison with what was done for Jacobi iterations; this is because of the overall greater relative
importance of SOR as a practical solution procedure. Second, it is written for general (nonsparse)
systems of linear equations. We have emphasized that it is almost always sparse systems that are
solved via iterative methods. For these, the algorithm simplifies somewhat. We will encounter cases
of this later in Chap. 5. Finally, we want to draw attention to the fact that there are really two
separate criteria by which the algorithm can be stopped: i) satisfaction of the iteration convergence
tolerance, and ii) exceeding the maximum permitted number of iterations, maxitr. The second
of these is crucial in a practical implementation because we do not know ahead of time whether
convergence to the required tolerance can be achieved. If it happens that it cannot be (e.g., due to
round-off errors), then iterations would continue forever unless they are stopped due to exceeding
the maximum specified allowable number.
1.2.3
In this subsection we will briefly summarize the foregoing discussions in the form of a table. presents
a listing of the main classes of problems one encounters in solving linear systems of equations. For
each type we give the preferred solution method, required storage for a typical implementation,
and total floating-point arithmetic needed to arrive at the solution, all in terms of the number
N of equations in the system. We note that the total arithmetic quoted for iterative solution of
sparse systems is a typical result occurring for SOR applied with optimal relaxation parameter
to solution of the discrete two-dimensional Poisson equation (see Chap. 5 for more details). In
25
general, for the case of sparse systems of this type, the arithmetic operation count can range from
slightly greater than O(N ) to as much as O(N 2 ), depending on the specific method employed.
A final remark on the comparison between direct and iterative methods, in general, is also in
order. Although the operation count is high for a direct method applied to a nonsparse system, this
count is precise: we can exactly count the total arithmetic, a priori, for such methods. Moreover,
this amount of arithmetic leads to the exact solution to the system of equations, to within the
precision of the machine arithmetic.
Table 1.1: Summary of methods for linear systems
By way of contrast, we can never exactly predict the total arithmetic for an iterative procedure
because we do not know ahead of time how many iterations will be required to achieve the specified
accuracy (although in simple cases this can be estimated quite well). Furthermore, the solution
obtained is, in any case, accurate only to the prescribed iteration tolerance, at bestdepending on
details of testing convergence. On the other hand, in many practical situations it makes little sense
to compute to machine precision because the problem data may be accurate to only a few significant
digits, and/or the equations, themselves, may be only approximations. All of these considerations
should ultimately be taken into account when selecting a method with which to solve a given linear
system.
1.3
The second main class of problems encountered in numerical linear algebra is the algebraic eigenvalue problem. As noted earlier, eigenvalue problems occur somewhat less frequently than does
the need to solve linear systems, so our treatment will be rather cursory. Nevertheless, in certain
areas eigenvalue problems are extremely important, e.g., in analysis of stability (in almost any
context) and in modal analysis of structures; hence, it is important to have some familiarity with
the treatment of eigenvalue problems. In this section we will begin with a brief description of the
eigenvalue problem, itself. We will then give a fairly complete treatment of the simplest approach
to finding eigenvalues and eigenvectors, the power method. Following this we briefly discuss what
is known as inverse iteration, an approach used primarily to find eigenvectors when eigenvalues are
already known. We then conclude the section with a short, mainly qualitative description of the
QR method, one of the most general and powerful of eigenvalue techniques.
Eigenvalue problems are of the form
AX = X ,
(1.34)
26
(1.35)
and only for such . It is not hard to check that (1.35) is a polynomial of degree N in if A
is a N N matrix. Thus, one method for finding eigenvalues of a matrix is to find the roots of
this characteristic polynomial. If N 4 this can be done exactly, although considerable algebraic
manipulation is involved for N > 2. Moreover, the eigenvectors must still be determined in order
to obtain a complete solution to the eigenvalue problem. We shall not consider such approaches
any further, and will instead direct our attention toward numerical procedures for calculating
approximate results for arbitrary finite N .
1.3.1
The power method is probably the simplest of all numerical methods for approximating eigenvalues.
As we will see, it consists of a fixed-point iteration for the eigenvector corresponding to the largest
(in magnitude) eigenvalue of the given matrix. Construction of the required fixed-point iteration
can be motivated in the following way.
Let {Xi }N
i=1 denote the set of eigenvectors of the N N matrix A. We will assume that these
eigenvectors are linearly independent, and thus form a basis for RN . Then any vector in RN can
be expressed as a linear combination of these eigenvectors; i.e., for any Y RN we have
Y = c1 X1 + c2 X2 + + cN XN ,
where the ci s are constants, not all of which can be zero. Multiplication of Y by the original matrix
A results in
AY = c1 AX1 + c2 AX2 + + cN AXN ,
(1.36)
But since the Xi s are eigenvectors, we have AXi = i Xi i = 1, 2, . . . , N . Thus (1.36) becomes
AY = c1 1 X1 + c2 2 X2 + + cN N XN .
Now if Y is close to an eigenvector of A, say Xi , then ci will be considerably larger in magnitude
than the remaining cj s. Furthermore, if we multiply by A a second time, we obtain
A2 Y = c1 21 X1 + c2 22 X2 + + cN 2N XN .
Clearly, if |1 | > |2 | > > |N | then as we continue to multiply by A, the term corresponding to
the largest eigenvalue will dominate the right-hand side; moreover, this process will be accelerated
if Y is a good estimate of the eigenvector corresponding to the dominant eigenvalue.
We now observe that the heuristic discussion just given does not supply a useful computational
procedure without significant modification because none of the Xi s or ci s are known a priori. In
fact, it is precisely the Xi corresponding to the largest that is to be determined. There are a
couple standard ways to proceed; here we will employ a construct that will be useful later. From
the original eigenvalue problem, AX = X, we obtain
X , AX = X T , X ,
27
X T , AX
.
=
hX T , Xi
(1.37)
The right-hand side of (1.37) is called the Rayleigh quotient. If X is an eigenvector of A, we can
use this formula to directly calculate the corresponding eigenvalue. We remark that (1.37) holds
for any operator A acting on elements of an inner-product space; i.e., if an operator and any one
of its eigenvectors (eigenfunctions in the infinite-dimensional case) are known, then the associated
eigenvalue can be immediately calculated.
We now employ (1.36) to construct a fixed-point iteration procedure to determine the eigenvector X corresponding to the largest eigenvalue, . We write
m
X (m) = AX (m1) Am Y = c1 m
1 X1 + + cN 1 XN
m
N
X
j
cj Xj .
= m
1
1
(1.38)
j=1
Now notice that in the Rayleigh quotient for finding , the true eigenvectors must be replaced by
the iterates defined by (1.38). Thus, the value obtained for is itself an iterate, defined by
E
D
T
X (m) , AX (m)
E .
(m+1) = D
T
X (m) , X (m)
But from (1.38) we have
AX (m) = X (m+1) .
Thus, we can write the iterate of simply as
D
(m+1)
E
T
X (m) , X (m+1)
E .
= D
T
X (m) , X (m)
(1.39)
We now want to show that this iteration procedure based on Eqs. (1.38) and (1.39) does, in
fact, converge to the largest eigenvalue of the matrix A. From (1.38) and (1.39) we have
(m+1)
N
X
T
i m j m+1
ci cj Xi , Xj
1
1
i,j=1
= 1 N
.
X m m
j
i
ci cj XiT , Xj
1
1
i,j=1
We have assumed from the start that |1 | > |2 | > > |N |, so as m the first term in the
series on the right will dominate. Thus it follows that
lim (m+1) = 1 + O ((2 /1 )m ) .
(1.40)
Hence, the sequence {(m) } converges to 1 as we expected. This also indirectly implies convergence
of the fixed-point iteration given by Eq. (1.38). In addition, it should be noted that in the case of
symmetric matrices, the order of the error term in (1.40) increases to 2m. (For details, see Isaacson
28
and Keller [15].) It is this feature that makes the present approach (use of the Rayleigh quotient)
preferable to the alternative that is usually given (see, e.g., Hornbeck [14]).
The preceding is a rather standard approach to arriving at the convergence rate of the power
method, but it does not motivate its construction as will appear in the algorithm given below. To
accomplish this we return to Eq. (1.34), the eigenvalue problem, itself, expressed as
X = AX .
We observe that this immediately takes on the fixed-point form if we simply divide by the eigenvalue
to obtain
1
X (m+1) = AX (m) .
We next observe that since eigenvectors are unique only to within multiplicative constants,
it is necessary to rescale the computed estimate after each iteration to avoid divergence of the
iterations; viz., lack of uniqueness ultimately results in continued growth in the magnitude of X if
normalization is not employed. Thus, we replace the above with
1
X (m)
X (m+1)
=
,
A
kX (m) k
kX (m+1) k
or
X (m+1) =
1 kX (m+1)k
AX (m) .
kX (m) k
(1.41)
This result appears rather different from the first equality in (1.38), but in fact, they are
equivalent. This follows from the fact that (up to an algebraic sign) it can be shown that the
eigenvalue is given by
kX (m+1) k
.
(1.42)
= lim
m kX (m) k
Indeed, this is often used in computational algorithms with k k taken to be the max norm because
this is far more efficient on a per iteration basis than is use of the Rayleigh quotient. On the other
hand, as alluded to above, use of the Rayleigh quotient for calculation of , as we have emphasized
herein, leads to higher convergence rates for the power method. In any case, use of (1.42) in (1.41)
leads to the usual fixed-point iteration form
X (m+1) = AX (m) ,
which when expanded as
X (m+1) = AX (m) = A2 X (m1) = = Am+1 X (0)
29
(m+1) =
5. If m = 0, go to 7
E
T
X (m) , X (m+1)
.
X (m)
2
2
1.3.2
The inverse iteration algorithm can be derived by viewing the algebraic eigenvalue problem, Eq.
(1.34), as a system of N equations in N + 1 unknowns. An additional equation can be obtained
by the typical normalization constraint needed to uniquely prescribe the eigenvectors, as discussed
above. This (N + 1) (N + 1) system can be efficiently solved using Newtons method (to be
presented in Chap. 2). Here we will merely present the pseudo-language algorithm for inverse
iteration, and refer the reader to Ruhe [25] for more details.
Algorithm 1.7 (Inverse Iteration with Rayleigh Quotient Shifts)
1. Set iteration tolerance and maximum number of iterations, maxitr;
Load initial guess for X and into X (0) , (0) .
Set iteration counter, m = 0.
30
5. Test convergence:
If
Y (m+1) Y (m)
+ |(m+1) (m) | < ,
then print results and stop.
1.3.3
The QR algorithm
The QR method is one of the most efficient, and widely used, numerical algorithms for finding all
eigenvalues of a general NN matrix A. It is constructed in a distinctly different manner from our
two previous methods. Namely, like a number of other procedures for computing eigenvalues, the
QR method employs similarity transformations to isolate the eigenvalues on (or near) the diagonal
of a transformed matrix. The basis for such an approach lies in the following concepts.
Definition 1.9 Let A and P be N N matrices with P being nonsingular. Then the matrix A
P AP 1 is said to be similar to the matrix A, and P ( )P 1 is called a similarity transformation.
A fundamental property of similarity transformations is contained in the following:
31
=
Theorem 1.5 The spectrum of a matrix is invariant under similarity transformations; i.e., (A)
N
(A), where denotes the set of eigenvalues {i }i=1 .
These ideas provide us with a very powerful tool for computing eigenvalues. Namely, we attempt
to construct a sequence of similarity transformations that lead to diagonalization of the original
matrix A. Since the new diagonal matrix has the same eigenvalues as the original matrix A (by
the preceding theorem), these can be read directly from the new (transformed) diagonal matrix.
The main difficulty with such an approach is that of constructing the similarity transformations.
There is nothing in the above definition or theorem to suggest how this might be done, and in
fact such transformations are not unique. Thus, the various algorithms utilizing this approach for
computing eigenvalues can be distinguished by the manner in which the similarity transformations
are obtained.
In the QR method this is done in the following way. It is assumed that at the mth iteration the
matrix A(m) can be decomposed as the product of a unitary matrix Q(m) and an upper triangular
matrix R(m) . (A unitary matrix is one whose inverse equals its transpose; i.e., Q1 = QT .) Hence,
A(m) = Q(m) R(m) .
Then we calculate A(m+1) as A(m+1) = R(m) Q(m) . It is easily checked that A(m+1) is similar to A(m)
and thus has the same eigenvalues. But it is not so easy to see that the eigenvalues can be more
easily extracted from A(m+1) than from the original matrix A. Because A need not be symmetric, we
are not guaranteed that it can be diagonalized. Hence, a fairly complicated procedure is necessary
in implementations of the QR method. The interested reader is referred to Wilkinson and Reinsch
[38] for details.
We have already noted that the QR method is very efficient, and that it is capable of finding
all eigenvalues of any given matrix. Its major disadvantage is that if eigenvectors are also needed,
a completely separate algorithm must be implemented for this purpose. (Of course, this can be
viewed as an advantage if eigenvectors are not needed!) Inverse iteration, as described in the
preceding section, is generally used to find eigenvectors when eigenvalues are computed with a QR
routine.
We will not include a pseudo-language algorithm for the QR method. As may be inferred from
the above discussions, such an algorithm would be quite complicated, and it is not recommended
that the typical user attempt to program this method. Well-validated, highly-efficient versions of
the QR procedure are included in the numerical linear algebra software of essentially all major
computing systems.
1.3.4
In the preceding subsections we have provided a brief account of methods for solving algebraic
eigenvalue problems. This treatment began with the power method, an approach that is straightforward to implement, and which can be applied to determine the largest (in magnitude) eigenvalue
of a given matrix and its corresponding eigenvector. The second method considered was inverse
iteration with Rayleigh quotient shifts. This technique is capable of finding an arbitrary eigenvalue
and its associated eigenvector provided a sufficiently accurate initial guess for the eigenvalue is
supplied. As we have already noted, this procedure is more often used to determined eigenvectors
once eigenvalues have already been found using some other method. One such method for finding
all eigenvalues (but no eigenvectors) of a given matrix is the QR method, the last topic of our discussions. This is a highly efficient technique for finding eigenvalues, and in conjunction with inverse
iteration provides a very effective tool for complete solution of the algebraic eigenvalue problem.
32
In closing this section we wish to emphasize that the methods described above, although widely
used, are not the only possible ones, and in some cases might not even be the best choice for a
given problem. Our treatment has been deliberately brief, and we strongly encourage the reader to
consult the references listed herein, as well as additional ones (e.g., the classic by Wilkinson [37],
and the EISPACK Users Manual [28]), for further information.
1.4
Summary
This chapter has been devoted to introducing the basics of numerical linear algebra with a distinct emphasis on solution of systems of linear equations, and a much briefer treatment of the
algebraic eigenvalue problem. In the case of solving linear systems we have presented only the
most basic, fundamental techniques: i) the direct methodsGaussian elimination for nonsparse
systems, tridiagonal LU decomposition for compactly-banded tridiagonal systems, and ii) the iterative schemes leading to successive overrelaxation, namely, Jacobi and GaussSeidel iteration, and
then SOR, itself. Discussions of the algebraic eigenvalue problem have been restricted to the power
method, inverse iteration with Rayleigh quotient shifts, and a very qualitative introduction to the
QR method.
Our principal goal in this chapter has been to provide the reader with sufficient information
on each algorithm (with the exception of the QR procedure) to permit writing a (fairly simple)
computer code, in any chosen language, that would be able to solve a large percentage of the
problems encountered from the applicable class of problems. At the same time, we hope these
lectures will provide some insight into the workings of available commercial software intended to
solve similar problems.
Finally, we must emphasize that these lectures, by design, have only scratched the surface of
the knowledge available for each chosen topicand barely even that for the algebraic eigenvalue
problem. In many cases, especially with regard to iterative methods for linear systems, there are far
more efficient and robust (but also much more elaborate and less easily understood) methods than
have been described herein. But it is hoped that the present lectures will provide the reader with
adequate foundation to permit her/his delving into the theory and application of more modern
methods such as Krylov subspace techniques, multigrid methods and domain decomposition. An
introduction to such material can be found in, e.g., Saad [26], Hackbusch [12], Smith et al. [27],
and elsewhere.
Chapter 2
2.1
In this section we will derive two fixed-point algorithms for single nonlinear equations: basic
fixed-point iteration and Newton iteration. The former is linearly convergent, in a sense to be
demonstrated below, while the latter is quadratically convergent.
2.1.1
Here we consider a single equation in a single unknown, that is, a mapping F : R R, and require
to find x R such that F (x ) = 0. We can, as in Chap. 1, express this in fixed-point form
x = f (x )
(2.1)
for some function f that is easily derivable from F . If we approximate f by a Taylor series expansion
about a current estimate of x , say x(m) , and retain only the first term, we have
x = f (x(m) ) + O(x x(m) ) .
Since the second term on the right cannot be evaluated without knowing x , we are forced to
neglect it and replace x on the left-hand side with a new estimate x(m+1) , resulting in the usual
fixed-point form
x(m+1) = f (x(m) ) ,
(2.2)
33
34
(2.3)
(2.4)
2.1.2
Newton iteration
Here, we begin by deriving the Newton iteration formula via a Taylor expansion, and we show that
it is in the form of a fixed-point iteration. Then we consider the convergence rate of this procedure
and demonstrate that it achieves quadratic convergence. Following this we view Newtons method
from a graphical standpoint, and as was the case in Chap. 1 with basic fixed-point methods, we will
see how Newtons method can failand we introduce a simple modification that can sometimes
prevent this. Finally, we present a pseudo-language algorithm embodying this important technique.
Derivation of the Newton iteration formula
We again consider the equation
F (x ) = 0 ,
(2.5)
F : R R. But now, instead of writing this in fixed-point form, we will employ a Taylor expansion
applied directly to the function F . For the mth approximation to x we have
F (x ) = F (x(m) ) + F (x(m) )(x x(m) ) + O((x x(m) )2 ) + .
If we now note that the left-hand side must be zero by Eq. (2.5), we obtain
F x(m)
(m)
x x
(m) + O((x x(m) )2 ) + .
F x
Applying arguments analogous to those used in the preceding section, we replace x with x(m+1)
on the left-hand side and drop the last term on the right-hand side to arrive at
F x(m)
(m+1)
(m)
(2.6)
x
=x
(m) .
F x
35
Clearly, if we define
G(x
(m)
)x
(m)
and write
F x(m)
(m)
F x
(2.7)
x(m+1) = G(x(m) ) ,
we see that Newtons method is also in fixed-point form; but the iteration function is more complicated than that for Picard iteration.
It can be shown that if F C 2 and F 6= 0 in an interval containing x , then for some > 0 and
x [x , x + ], the Lipschitz constant for G is less than unity, and the sequence generated by
(2.6) is guaranteed to converge to x . In particular, we have
F x(m) F x(m)
F (x(m) )
(m)
G (x ) = 1 (m) +
.
2
F x
F x(m)
(m+1)
(m)
x x
=x x
+ (m) ,
F x
1
0 = F (x ) = F (x(m) ) + F (x(m) )(x x(m) ) + F (x(m) )(x x(m) )2 + .
2
Thus,
1
F (x(m) ) = F (x(m) )(x x(m) ) F (x(m) )(x x(m) )2 .
2
We substitute this into the above, and use the definition of em , Eq. (2.4), to obtain
1 F x(m) 2
e .
em+1 = em em
2 F x(m) m
Hence, with F C 2 and F (x(m) ) 6= 0, this reduces to
|em+1 | = K|em |2 ,
(2.8)
36
with K < . So the error at the (m + 1)th iteration of Newtons method is proportional to the
square of the error at the mth iteration. This gives rise to the often used terminology quadratic
convergence of Newtons method. It is easy to see from (2.8) that once the error is O(1) it decays
very rapidly. In fact, once the error is reduced to single precision machine , only one additional
iteration is required to achieve double precision accuracy. This very rapid convergence is one of the
main reasons for the wide use of Newtons method.
Graphical interpretation of Newtons method
It is worthwhile to examine Newtons method from a geometric standpoint, just as we did for fixedpoint iteration in Chap. 1. Consider the following figure depicting a sequence of Newton iterations.
Geometrically, the iterations proceed as follows.
Starting from an initial guess x(0) , evaluate
F at
(0)
(0)
(0)
(0)
x . Then from the point x , F (x ) draw a tangent corresponding to F x
to intersect
the x-axis. This point of intersection is the next iterate, x(1) . It is of interest to note that if F is
linear, Newtons method converges ina single iteration, independent of the initial guess, since the
slope through the point x(0) , F (x(0) ) is the graph of F , itself. Furthermore, it is easily seen from
the graphical construction that convergence is very rapid, even when F is not linear, in agreement
with what we would expect from (2.8).
F(x)
F(x (0) )
x (1)
F(x (2) )
x (2)
x*
x (0)
F(x (1) )
37
rapidly. This unfortunate behavior occurs in this particular problem simply because the Newton
step,
F (x)
,
x =
F (x)
is too large, since F is very small and F 1. The remedy is to replace x with x, with
0 < 1.
F(x)
F(x (0) )
x (1)
x*
x (0)
F(x (1) )
F x(m)
x = (m) ,
F x
x(m+1) = x(m) + x .
38
We note that convergence tests other than the one employed in the above algorithm may be used.
Philosophically, at least, we would expect to have |F (x(m) )| < since we are seeking x F (x ) = 0.
Nevertheless, the test given above is probably the most widely used because one can argue that
once x becomes sufficiently small, additional corrections will have no effect, and will thus be a
waste of arithmetic. But the following should be noted. Clearly, if |F | is very small but nonzero
close to the desired root, the test used in the above algorithm may be unreasonably conservative.
On the other hand, if |F | |F |, convergence may be indicated for an inaccurate approximation
to the root. It is thus sometimes useful to employ a test of the form max (|x|, |F |) < .
2.2
One of the main objections to the use of Newtons method has always been the requirement to
calculate the derivative of the function F whose zero is sought. In some situations this is merely
tedious, but in others either the derivative does not exist at all, or it cannot be expressed in terms
of elementary functions. In the former case, it may be advisable to reformulate the problem so
that the new function is differentiable, or possibly an iteration scheme not requiring derivatives (or
their approximations) might be used. On the other hand, if difficulty in evaluating the derivative
of F is the only problem, the methods to be presented in this section should be of value.
2.2.1
These first of these methods to be considered is the secant method. To construct the iteration
function for this method we begin with the Newton iteration formula, Eq. (2.6):
F x(m)
(m+1)
(m)
x
=x
(m) ,
F x
and replace F with an approximation that depends only on values of F , itself. In the next chapter we will consider several such approximations in fair detail. The secant method is typically
constructed using a simple, intuitive one. In particular, we use
F x(m) F x(m1)
(m)
+ O(xm x(m1) ) ,
F (x ) =
x(m) x(m1)
and the iteration formula becomes
x
(m+1)
=x
(m)
x(m) x(m1) F x(m)
F x(m) F x(m1)
(2.10)
= x(m) + x .
There are several things to notice with regard to this iteration scheme. The first is that it is
not a fixed-point method because the iteration function depends on more than just the most recent
39
F(x)
F(x (1) )
F(x (2) )
x (0)
x (3)
x* x (4)
x (2)
x (1)
F(x (3) )
F(x (0) )
many of the difficulties encountered with Newtons method can also occur with the secant method,
and that one of the remedies is use of a damping factor. In particular, we can replace x with
x in Eq. (2.10), as was done in Eq. (2.6) for Newtons method. An additional problem with the
secant method, not present in Newtons method, occurs as the zero of F is approached. Namely,
F (x(m) ) F (x(m1) ) can approach zero, leading to floating-point overflow during evaluation of Eq.
(2.10). The main remedy for this problem is use of higher-precision machine arithmetic.
The algorithm for the secant method is very similar to the Newton algorithm, with the exception
of the definition of x. Hence, we shall not present a formal pseudo-language treatment for this
case, and instead leave this construction as an exercise for the reader.
2.2.2
For one reason or another, it is sometimes necessary to bound the interval on which a root of an
equation is sought. Usually, but not always, this is because the equation has more than one zero,
and a particular one may be desired. For example, if we wish to find all zeros of a polynomial of
degree n and we have already found some of them, we would want to restrict the intervals on which
40
further roots are sought. One procedure for doing this is the method of false position, or regula
falsi. This method is rather similar to the secant method in that it does not require evaluation
of derivatives. But at the same time, in contrast to the secant method, the derivative appearing
F(x)
F(x (1) )
F(x (2))
(3)
x (0) F(x )
x* x (3)
x (2)
x (1)
F(x (0) )
(b a)F (b)
.
F (b) F (a)
41
5. Calculate x(m+1) :
If F (x(m) ) F (a) < 0, then
(m+1)
=x
If F x(m) F (b) < 0, then
x
6. If m < maxitr, then m = m + 1
(m+1)
(m)
x(m) a F x(m)
.
F x(m) F (a)
b x(m) F (b)
.
=b
F (b) F x(m)
go to 4
else print error message and stop
An interesting point to make regarding this algorithm is that, unlike the secant method, this is
a fixed-point iteration with iteration function
x(m) a F x(m)
(m)
if F x(m) F (a) < 0 ,
x
(m)
F x
F (a)
(m)
(2.11)
=
f x
(m)
bx
F (b)
if F x(m) F (b) < 0 .
b
F (b) F x(m)
As might be expected from earlier discussions, it can be shown that this iteration function results
in a linear convergence rate. Hence, it is not very efficient. On the other hand, if a zero exists on
a specified interval, and F is continuous, the regula falsi method is guaranteed to find it.
2.3
In Chap. 1 we derived a basic fixed-point iteration scheme to treat systems of equations from
the start because we wished to apply it to linear systems. In the present chapter we began by
considering only a single nonlinear equation because convergence rate estimates are more easily
obtained for this case. But in practice we very often must face the task of solving systems of
nonlinear equations. Rather typically, Picard iteration is not suitable because either convergence
is too slow, or it does not occur at all. Because of this, Newtons method is the most widely-used
approach for treating nonlinear systems. It can be shown (via the NewtonKantorovich theorem,
see e.g., Stoer and Bulirsch [30]) that Newton iterations still converge quadratically for nonlinear
systems, just as for a single equation. Here, we will present a detailed construction of the method
and a pseudo-language algorithm for its implementation.
2.3.1
While this method can be derived via a formal Taylor series expansion as was used to obtain the
1-D iteration earlier, in multi-dimensional cases this is somewhat tedious, and we will here simply
argue from analogy with the 1-D iteration formula, Eq. (2.6).
42
provided F 6= 0 and F is bounded for x near x . To generalize this to the case N > 1, we need
a reasonable interpretation of the derivative of mappings F : RN RN . It turns out that the
Jacobian matrix of F is the correct interpretation of the derivative for this case. Hence, if we let
x = (x1 , x2 , . . . , xN )T and F = (F1 , F2 , . . . , FN )T , then the Newton iteration scheme becomes
(m+1) (m)
F1 1 (m)
F1 F1
F1
x1
x1
x1
x2
xN
F2 F2
F2
F2
x2
x2
x
x2
xN
1
(2.12)
=
FN FN
FN
FN
xN
xN
x1
x2
xN
In more concise notation this is
i1
h
F (m)
x(m+1) = x(m) J F (m)
(2.13)
= x(m) + x ,
where J(F ) denotes the Jacobian matrix of F , shown in detail on the right-hand side of (2.12).
From this it follows (by definition) that
h
i1
x J F (m)
F (m) .
It should be noted that, contrary to what might at first appear to be the case, it is not necessary to
compute the inverse Jacobian matrix
when using Newtons method. In particular, upon multiplying
the above expression by J F (m) we obtain
(2.14)
J F (m) x = F (m) .
Once x(m) is knownfor example, x(0) is specified at the startF (m) can be evaluated (this
is just the system to be solved), and J F (m) can be directly calculated since this simply consists
of the partial derivatives of F evaluated at x(m) . Thus, Eq. (2.14) is a linear system for x, with
known system matrix and right-hand side. It can be solved using Gaussian elimination described
in Chap. 1. After x has been determined, we can immediately obtain the updated approximation
x(m+1) to x from the second line in (2.13). We see from this that solving systems of nonlinear
equations mainly requires that we only be able to solve systems of linear equations.
We note that, in general, damping may be required just as for a single equation. But it can now
be introduced in several ways. First, one can use the same amount of damping for every equation
by multiplying each component of x by the scalar constant 0 < 1. One might also employ
a damping vector (formally, a diagonal matrix) which, in general, would allow damping of each
component of the solution by a different amount. Finally, even a damping matrix, say D, might be
utilized. Although the latter two alternatives are seldom used, it is important to be aware of these
possibilities because they can be valuable in special situations.
43
2.4. SUMMARY
2.3.2
We close this section with the pseudo-language algorithm for applying Newtons method to systems
of nonlinear equations.
Algorithm 2.3 (Newtons Method for Systems)
1. Input , and maxitr; initialize solution vector x(0) ; set m = 0.
2. Evaluate the vector F (x(m) ).
3. Test convergence:
If kF k < , then print results and stop.
else evaluate elements of Jacobian matrix, J F (x(m) ) .
J F (m) x = F (m)
2.4
Summary
In this chapter we have considered several of the main methods in wide use for solving nonlinear
algebraic equations and systems. The procedures treated here included basic fixed-point iteration
already encountered in Chap. 1 for linear equations, Newtons method, the secant method and
regula falsi for single equations, and Newtons method for systems. It is clear from the above
discussions that the basic fixed-point iterations represent the simplest approach, but they are not
usually preferred because of their rather poor convergence rate. Newtons method is probably the
most popular technique for solving nonlinear equations because of its quadratic convergence rate.
But it sometimes must be damped if poor initial guesses are used. It should also be mentioned
that Newtons method is sometimes started with a Picard iteration to improve the initial guess.
The secant method is formally the most effective of the approaches we have considered here in the
sense that its convergence rate is nearly that of Newtons method, but it requires only a single
44
function evaluation per iteration. Regula falsi provides the advantage of bounding the domain in
which a solution is sought, but as noted above, it is a basic fixed-point iteration; so it converges
only linearly.
We have treated only Newtons method for systems of nonlinear algebraic equations. It was
seen that this is a natural extension of the single equation case, and this along with its quadratic
convergence rate contribute to its popularity. Its main disadvantage is its requirement of calculating Jacobian matrix elements. There are a number of so-called quasi-Newton methods (which,
in a sense, extend the secant method to systems) that are often used, especially in the context of
nonlinear optimization, to alleviate this difficulty. However, it is seldom that full quadratic convergence rates are achieved by these techniques. We refer the interested reader to [30] for further
discussions on this topic.
Finally, we wish to emphasize that there are numerous other procedures for solving nonlinear
algebraic equations, and the choice of methods presented here at least in part reflects a bias of
the author. Conspicuous by its absence is the bisection method. This is usually the first method
treated in elementary texts because of its intuitive appeal. But its convergence rate is very low,
and it is quite difficult to extend to use for systems of equations; hence, we have chosen to ignore it
in these lectures. We have also neglected to discuss methods designed specifically for finding zeros
of polynomials. The methods we have presented are capable of doing this, but other approaches
can be far more efficient in specific situations. The reader should consult books devoted specifically
to solution of nonlinear equations by, for example Ostrowski [22] or Traub [36] for information on
various of these important and interesting subjects.
Chapter 3
Approximation Theory
In numerical analysis there are three main mathematical objects which we must often approximate:
i) functions, ii) integrals of functions, and iii) derivatives of functions. We have already seen a case
of the last of these in contructing the secant method for solving nonlinear equations. In this chapter
we will treat all three in some detail; but in all cases we will study only rather basic and widely used
methods. In the approximation of functions we consider least-squares methods briefly, and then
devote some time to exact Lagrange interpolation and construction of cubic splines. Approximation
of integrals will be done via the trapezoidal rule, Simpsons rule and GaussLegendre quadrature
for integrals on subsets of the real line. We then show how these can be applied recursively to
evaluate multiple integrals. Following this we will present methods for approximating derivatives,
mainly via the finite-difference approach. In a final section we consider one of the most important,
yet too often neglected, topics in numerical analysis, namely grid function convergence. We will
see that, at least in a limited way, we are able to assess the accuracy of our approximations based
only on the computed results, themselves.
3.1
Approximation of Functions
We begin this section with the least-squares method, for the most part because this provides a
good application of Newtons method for systems of equations, which was just presented at the end
of the preceding chapter. Then Lagrange interpolation is discussed. This method is widely used,
especially when abscissas are not equally spaced. We then discuss cubic spline interpolation. This
will provide the first concrete application of tridiagonal LU decomposition presented in Chap. 1.
3.1.1
Least-squares methods are among the most widely-used techniques in the analysis of experimental
data. In later sections we will consider the so-called exact approximation methods which force
the interpolation function to exactly equal the data at certain prescribed points. This is generally
inappropriate in analyses of data from almost any experiment since such data are usually of only
limited accuracy. Least-squares methods will not generally lead to exact reproduction of the data
value at any particular point. On the other hand, they provide the best possible fit of all the data
in a certain global sense.
The reader may have previously encountered the use of least-squares procedures under the name
linear regression, or multiple regression. The linear regression formulas are derived from the same
basic principles we will use here. But when data are to be fit to linear functions, specific use can be
45
46
made of linearity to simplify the computational formulas. In the present treatment we will assume
the data are to be fit to nonlinear functions (i.e., nonlinear in the parameters to be determined),
and no simplifications will be made. Our method will still be valid in the linear case, but it may
be computationally less efficient than would be use of linear regression formulas.
Some General Practical Considerations
As a starting point we consider the graph of some experimental data, displayed in Fig. 3.1. It
may be desirable for any of a number of reasons to find a single function with which to represent
these data. The first step in the general nonlinear case is to select a functional form that exhibits
many of the qualitative features of the data. It is always a good idea to first plot the data to gain
some insight into the choice of function for the attempt at curve fitting. In the present case, the
data are somewhat suggestive of a damped sine wave, but with phase shift, and possibly x-varying
frequency. We could no doubt fit these data reasonably well by utilizing a polynomial of degree at
least three (why degree three?). Often there may be theoretical reasons for choosing a particular
function with which to fit the data; but in the absence of any theory, a polynomial of sufficient
degree to accommodate all of the zeros and local maxima and/or minima is a reasonable choice.
For the data in Fig. 3.1, we might on intuitive grounds choose
f ( x)
Experimental data
Leastsquares fit
(3.1)
to perform the least-squares correlation. In this expression , , , A and B are all unknown
parameters to be determined from the data via the least-squares method. It is not unusual to fit
data to this number of parameters, and it is worthwhile in the present case to examine just what
influence each of these has on the value of the function g. The parameter A is the amplitude of
the sine function with respect to the mean value, B. For the data shown in Fig. 3.1 we would
probably estimate B 0 and A less than or equal to the absolute maximum data point value. The
parameter is the exponential decay rate, which can be easily estimated from the data; is the
phase shift, which should be approximately /2 in the present case. Finally, is the frequency,
which again can be estimated from specific data. It is important to be able to make reasonably
47
good estimates of the parameters for a nonlinear curve fit of data because these are needed as
initial guesses in the iteration scheme (usually Newtons method) that must be employed to solve
the nonlinear equations whose solutions are the parameter values.
Least-Squares Method Construction
We now consider more details of how to determine the parameters of the least-squares fit. We
first need to make clear what it is that we wish to accomplish. Suppose there are N data pairs
(xi , f (xi )), i = 1, . . . , N that we wish to represent by the function g. This function, of course,
depends on x, but it also depends on the five unknown parameters , , , A and B. Let the
vector = (1 , 2 , . . . , M )T formally represent these parameters, and express g as g(x, ). The
goal of a least-squares correlation (fit) of data is to select the components i of the parameter vector
such that the sum of the squared differences [g (xi , ) f (xi )]2 over all data points is as small as
possible. In other words, we want to choose to minimize
S(x, ) =
N
X
i=1
(3.2)
S(x, ) is often referred to as the least-squares functional, and we should observe that it is actually
a function of only the parameter vector that is, dependence on x has been summed out.
We know from vector calculus that the minimum of a function occurs at the point where its
gradient is zero, and at the same time its Hessian matrix (the matrix of second partial derivatives of
S()) is positive definite. We shall assume, for convenience that the latter holds, but we note that
this is not automatic. We now view S as a function only of since, as already noted, performing
the summation removes the x dependence. Then the minima of S occur for those such that
S = 0 ;
that is for such that
S
= 0,
j
j = 1, 2, . . . , M .
(3.3)
Equations (3.3) comprise a system of M nonlinear equations called the normal equations, for the
M (< N ) unknown parameters, i (in the case of Eq. (3.1), M = 5). Thus, we expect to be able to
solve these via Newtons method. Define
S
, j = 1, 2, . . . , M ;
Fj (1 , 2 , . . . , M )
j
then for F = (F1 , . . . , FM )T we have
i1
h
(m+1) = (m) J F (m)
F (m) .
(3.4)
X g
S
=2
[g (xi , ) f (xi )] = 0,
j
j
i=1
j = 1, 2, . . . , M .
(3.5)
A second differentiation of this result leads to the Jacobian matrix required for application of
Newtons method, which is the Hessian matrix of the least-squares functional:
N
X
2g
2S
g g
Fj
=
=2
(g(xi , ) fi ) +
,
j, k = 1, 2, . . . , M .
(3.6)
k
j k
j k
j k
i=1
48
Detailed Example
Our main task in constructing the least-squares formulation becomes one of finding the analytical
expressions for the components of F and J(F ). We will carry out part of this calculation for the
function given in Eq. (3.1) to provide a concrete demonstration of the technique.
Now in our example function given in Eq. (3.3), 1 = , 2 = , 3 = , 4 = A, and 5 = B.
Thus,
g
g
=
1
g
g
=
2
g
g
=
3
g
g
=
4
A
g
g
=
5
B
We now use the common shorthand notation f (xi ) = fi that is used extensively in numerical
analysis and write the system of five equations to be solved for = (, , , A, B)T by Newtons
method:
F1 () =
F2 () =
N n
X
i=1
N n
X
i=1
F3 () =
F4 () =
F5 () =
N n
X
i=1
N n
X
i=1
N n
X
i=1
io
h
xi Aexi cos (xi + ) Aexi sin (xi + ) + B fi = 0
io
h
xi Aexi sin (xi + ) Aexi sin (xi + ) + B fi = 0
io
h
Aexi cos (xi + ) Aexi sin (xi + ) + B fi = 0
io
h
exi sin (xi + ) Aexi sin (xi + ) + B fi = 0
o
Aexi sin (xi + ) + B fi = 0 .
The next step is to derive the Jacobian matrix J(F ) of the system, corresponding to Eq.
(3.6). We already have the expressions required to construct the second term in brackets in these
equations. All that remains is to derive the second partial derivatives of g. We will carry this out
49
2g
2g
A
2g
B
We leave the remaining calculations to the reader, and note that this might typically be done using
symbolic manipulation languages such as Maple or Macsyma, for example.
The algorithm for solving least-squares problems merely consists of the Newtons method algorithm of Chap. 2 in which the nonlinear equations are evaluated using Eq. (3.5) (with the factor of
2 deleted), and Eqs. (3.6) are employed to construct the required Jacobian matrix. Thus, we shall
not present additional details; but it is worth mentioning that the same procedure can be employed
for minimization of general nonlinear functions. Its application is not restricted to least-squares
problems.
3.1.2
The next method to be considered for function approximation is Lagrange interpolation. We should
all be familiar with linear interpolation; we will see that this is the simplest case of Lagrange interpolation, one of a class of methods termed exact. This terminology comes from the fact that
Lagrange polynomials are constructed so that they exactly reproduce the function values at the
abscissas of the given data points. For example, recall that the linear interpolation formula is
constructed so that it exactly reproduces both endpoint function values on the interval of interpolation because it is simply a straight line through these two points. Hence, it is a first degree exact
interpolation polynomial, and we will see below that it is a Lagrange polynomial.
We begin this section with an heuristic approach to exact interpolation that can sometimes
be useful in situations more general than the polynomial case to be treated, and we follow this
with a theorem (restricted to the polynomial case)on uniqueness of such results. We then present
construction of the Lagrange interpolation polynomials, analyze their accuracy and extend the
constuction procedure to two independent variables. Finally, we briefly remark on some weaknesses
of this form of interpolation.
A Basic Exact Interpolation Construction Method
It is probably of some value to build an exact interpolation polynomial in an intuitive (although
computationally inefficient) manner first in order to demonstrate some of the basic ideas. Suppose
we are given three points (x1 , f1 ), (x2 , f2 ) and (x3 , f3 ), and we are required to construct an exact
interpolation polynomial pn (x) of degree n which takes on the values f1 , f2 , f3 at x1 , x2 , x3 ,
respectively. Our first task is to determine the appropriate degree, n, of the polynomial with which
we intend to represent these data. In general, a polynomial of degree n is of the form
pn (x) = a0 + a1 x + a2 x2 + + an xn .
(3.7)
50
Thus, it contains n + 1 coefficients to be determined, and this implies that n + 1 points (xi , fi ) will
be required. Conversely, if n + 1 points are given, and it is required that they be interpolated with
an exact polynomial, the polynomial must in general be of degree at least n. In our particular case
we are given three points, so we must employ a polynomial of degree at least two, i.e., a quadratic.
The next step is to formally evaluate Eq. (3.7) at each of the data points. We obtain
p2 (x1 ) = f1 = a0 + a1 x1 + a2 x21
p2 (x2 ) = f2 = a0 + a1 x2 + a2 x22
(3.8)
p2 (x3 ) = f3 = a0 + a1 x3 + a2 x23
with the left-hand equalities arising due to the requirement of exactness. Now since the fi s and xi s
are given, we see that (3.8) is merely a linear system of three equations for the three ai s. It can be
demonstrated that the system matrix (sometimes called the Vandemonde matrix) is nonsingular if
and only if xi 6= xj when i 6= j ; i.e., the data are not multivalued. Thus, the solution exists and is
unique, and this implies the following:
1
Theorem 3.1 Let {xi }n+1
i=1 be a set of distinct points on R , and let fi = f (xi ). Then a unique
1
1
polynomial pn : R R of degree n,
pn (x) =
n
X
ai xi ,
i=0
such that
pn (xi ) = fi
i = 1, 2, . . . , n + 1 .
(3.9)
Several remarks are of value at this time. The first is that polynomials of higher degree can
be constructed through the same set of points merely by deleting one or more, as required, lowerdegree terms and replacing these with higher-degree ones. Second, although the system matrix
correponding to (3.8) can be shown to be nonsingular, it is often ill-conditioned. So coefficients
obtained as solutions to (3.8) may be subject to large rounding errors. Finally, although we will
restrict our attention to exact polynomial fits in these lectures, the procedure just described can
be used with other functional forms, e.g., rational functions, exponentials, etc.
Construction of Lagrange interpolation polynomials
We will now present a construction method that is numerically stable, and more efficient than
Gaussian elimination for the present application. To begin we assume that pn (x) can be expressed
as
n+1
X
i (x)fi ,
(3.10)
pn (x) =
i=1
where the fi s are the given ordinates. The problem now reduces to determining the i s. Clearly,
at least one of these must be a polynomial of degree n (and none can have degree > n) if this is to
be true of pn , and in Lagrange interpolation this is true for all the i s. Moreover, if pn (xi ) = fi is
to hold, we should expect that
(
1 for i = j
(3.11)
i (xj ) =
0 for i 6= j;
or
i (xj ) = ij ,
51
n+1
Y
j=1
j6=i
(x xj )
,
(xi xj )
i = 1, 2, . . . , n + 1.
(3.12)
It is easily checked that when x = xi , each factor in the product is unity, so i (xi ) = 1. Similarly
for x = xj for some j 6= i, one of the factors appearing in Eq. (3.12) will be zero, and we have
i (xj ) = 0. Finally, if x 6= xi i = 1, . . . , n + 1, all of the factors are nonzero. There are n such
factors, because the one corresponding to j = i is deleted. Hence, each i is of degree n.
The pseudo-language algorithm for Lagrange interpolation is easily constructed (simply evaluate
Eq. (3.12) for any desired value of x, and substitute the results into Eq. (3.10)), so we leave this to
the reader. However, the following remarks are of value. Each i requires 2n1 multiplies (actually,
n divisions and n 1 multiplies) for its construction. Then an additional n + 1 multiplies and n
adds are needed to evaluate Eq. (3.10). In any case the total arithmetic required to evaluate pn at
any desired point is only O(n2 ) by this method. If we were to follow the intuitive idea of solving
Eq. (3.8) for the polynomial coefficients ai by Gaussian elimination, O(n3 ) arithmetic operations
would be required just to determine the ai s. Additional work would still be needed to evaluate
pn (x). However, it should be mentioned that in many applications n is quite small; in these cases
there may be situations in which there is no distinct advantage in employing (3.10) and (3.12)
beyond the fact that they are general formulas that can be used for all n < .
A specific examplelinear interpolation
We now present a specific example of applying the preceding Lagrange interpolation formulas. To
provide a familiar context, and at the same time avoid undue algebraic complications, we will
construct the formula for linear interpolation between the two points (x1 , f1 ) and (x2 , f2 ). Since
there are two points, the Lagrange polynomial must be of degree one; i.e., it is a (generally affine)
linear function.
From (3.10), we have for n = 1
p1 (x) =
2
X
i=1
n+1
Y
j=1
j6=i
2
Y
(x xj )
(x xj )
x x2
1 (x) =
=
,
(xi xj )
(x1 xj )
x1 x2
j=1
j6=1
and
2 (x) =
2
Y
(x xj )
x x1
=
.
(x2 xj )
x2 x1
j=1
j6=2
Hence,
p1 (x) =
x x1
x x2
f1 +
f2 ,
x1 x2
x2 x1
(3.13)
52
f2 f1
(x x1 ) .
x2 x1
(3.14)
The form of (3.14) is equivalent to the standard one which would be obtained by writing the
equation of the straight line between (x1 , f1 ) and (x2 , f2 )and also what would be found via
Taylor expansion of f about x = x1 followed by replacing the derivative with a first divided
difference, hence, a Taylor polynomialthus suggesting uniqueness of at least the first-degree
Lagrange interpolation polynomial.
Accuracy of Lagrange polynomial interpolation
Another thing that can be deduced from Eq. (3.14) is the accuracy achieved by linear interpolation. Suppose we were given only the point (x1 , f1 ), but in addition were also given the first two
derivatives of f evaluated at x1 , i.e., f (x1 ) and f (x1 ). Then for any x close to x1 we can expand
f in a Taylor series about x1 :
1
f (x) = f (x1 ) + f (x1 )(x x1 ) + f (x1 )(x x1 )2 + .
2
(3.15)
We notice that this is the same as (3.14) through the first two terms, except that f (x1 ) in (3.15)
is replaced by the slope (f2 f1 )/(x2 x1 ) in (3.14). Now if we approximate f (x1 ) by a forward
difference (similar to what was used to develop the secant method in Chap. 2), we obtain
f (x1 ) f1 =
f2 f1
1
f1 (x2 x1 ) + ,
x2 x1 2
as we will show in more detail later in this chapter. Substitution of this into Eq. (3.14) yields
f (x) = f1 +
1
f2 f1
(x x1 ) + f1 (x x1 )2 (x2 x1 )(x x1 ) + .
x2 x1
2
(3.16)
We see that f (x) in (3.16), which is a precise (infinite Taylor expansion) representation, agrees
with p1 (x) up to the third term. This term is called the dominant truncation error, and since
x x2 , the term is bounded by C(x2 x1 )2 , where C is a constant. That is,
f (x) p1 (x) = C(x2 x1 )2 .
From this we see that linear interpolation is second-order accurate in the step size x = x2 x1 .
Thus, as x is reduced p1 (x) more accurately approximates f (x) x [x1 , x2 ], in agreement with
our intuition. But beyond this we see that the error in this approximation decreases with the square
of x.
Similar analyses show that, in general, a Lagrange interpolation polynomial of degree n provides
an approximation of order n + 1. In particular, it is shown (for example, in Isaacson and Keller
[15]) that if f (x) is approximated by a polynomial of degree n, then
!
n+1
Y
1
(x xi ) f (n+1) () ,
f (x) pn (x) =
(n + 1)!
i=1
where the xi are the abscissas of the data points used to construct pn , f (n+1) denotes the (n + 1)th
derivative of f , and
min(x1 , . . . , xn+1 , x) < < max(x1 , . . . , xn+1 , x) .
53
we have
f (x) pn (x) Cxn+1 O(xn+1 ) .
Thus, in general, a nth -degree Lagrange polynomial provides a (n + 1)th -order approximation.
Linear interpolation in two dimensions
We conclude this section on Lagrange interpolation by demonstrating how to apply linear interpolation in two independent variables. Rather than present a theoretical development, we will employ
an intuitive approach. In Fig. 3.2 we depict a grid of function values corresponding to a function
f (x, y): R2 R1 . Suppose we are given the values of f at the four points (x1 , y1 ), (x2 , y1 ), (x1 , y2 )
and (x2 , y2 ), and we are required to find a linear approximation to f at the point (x , y ). The
y2
y*
Interpolate f 1 and f 2
*
*
to obtain f**
y1
x* x2
f21 f11
(x x1 )
x2 x1
54
f22 f12
(x x1 )
x2 x1
f2 f1
(y y1 )
y2 y1
As already indicated, these expressions can be combined, analytically, to obtain the above
mentioned bilinear interpolation formula:
f = f11 +
f12 f11
f22 f21 f12 + f11
f21 f11
(x x1 ) +
(y y1 ) +
(x x1 ) (y y1 ) .
x2 x 1
y2 y1
(x2 x1 )(y2 y1 )
We leave as an exercise to the reader showing that this expression provides a second-order approximation to f (x, y) in the same sense already discussed in the one-dimensional case. (Hint: note the
resemblance of the preceding expression to a multi-dimensional Taylor expansion of f .) Finally,
we mention that even higher-dimensional linear approximations can be constructed in a manner
analogous to that summarized in Algorithm 3.1, and furthermore, higher-degree interpolation polynomials might be used in place of the linear ones treated here.
Difficulties with high-degree interpolation polynomials
We conclude this section by observing that one of the major shortcomings of using Lagrange
polynomials is that the high-degree polynomials required for formal accuracy often exhibit very poor
behavior away from the abscissas, as is illustrated in Fig. 3.3. Over long intervals it is impossible
pn ( x)
f (x)
x
Figure 3.3: Ill-behavior of high-order Lagrange polynomials
to obtain reasonable accuracy with a single low-degree polynomial; but the wiggles that occur
when high-degree global approximations are used may also be unacceptable. For example, it is
clear from Fig. 3.3 that pn (x) is not generally a good approximation to f (x)at some points even
its sign is incorrect, and this often can be a severe disadvantage. In the next section we will present
a partial remedy.
55
3.1.3
We probably know from experience that linear interpolation usually works reasonably well. The
reason for this it is that is usually applied over short subintervals of the domain of definition of
the function being interpolated. This is often sufficient when only function values are needed. But
if derivatives of the interpolated function are required, for example to be employed in Newtons
method, local linear interpolation should not be used. This can be seen from Fig. 3.4. In particular,
the derivative of the interpolation polynomial p1 (x) for x [x1 , x2 ] is exactly the constant (f2
f1 )/(x2 x1 ). On the other hand for x [x2 , x3 ] the derivative is (f3 f2 )/(x3 x2 ). Hence, there
is a discontinuity in p1 at x = x2 .
f (x )
p1(x) ,
x [ x2 ,x3 ]
p1(x) ,
x1
x2
x [ x1 ,x2 ]
x3
56
The first step in deriving the equations for determining the coefficients for spline interpolation
is to recognize that since the spline is a cubic on each subinterval, its second derivatives must be
x xi+1
x xi
Si+1
Si ,
hi+1
hi+1
where hi+1 xi+1 xi , and x [xi , xi+1 ]. Observe that by definition S is required to be continuous
x [a, b]. Note also that we do not assume the hi s to be equal. Next, we integrate the above
equation twice: the first integration yields
(x xi+1 )2
(x xi )2
Si+1
Si + Ai ,
2hi+1
2hi+1
(3.17)
(x xi+1 )3
(x xi )3
Si+1
Si + Ai (x xi ) + Bi .
6hi+1
6hi+1
(3.18)
S (x) =
and a second integration leads to
S(x) =
(xi xi+1 )3
S i + Bi = f i
6hi+1
1 2
h S + Bi = fi .
6 i+1 i
+ hi+1 Ai + Bi = fi+1 .
S(xi+1 ) = h2i+1 Si+1
6
Hence, the integration constants are
fi+1 fi hi+1
(Si+1 Si ) ,
hi+1
6
h2
Bi = fi i+1 (Si ) .
6
Ai =
(3.19)
If we substitute these into (3.18) and use the fact that (x xi+1 )3 = (x xi hi+1 )3 , after some
algebra we arrive at a canonical expression for S that is analogous to p1 (x) in (3.14):
fi+1 fi 1
S(x) = fi +
57
(3.20) is in the form of a truncated Taylor series for f (x), expanded about xi . From earlier results
we would expect
f (x) S(x) = O(h4 )
to hold, since S is locally a cubic. It is clear from (3.20) that in order for this to be true, the
following must be satisfied:
fi+1 fi 1
fi
= O(h) .
hi+1
(3.22)
(3.23)
It is easily shown that (3.22) implies (3.21) and (3.23); we leave this as an excercise for the reader.
However, (3.22) is not automatic. The reader is referred to deBoor [5] for details of this, and in
general, for extensive treatment of cubic splines.
Computational Algorithm
We now turn to development of a computational scheme for determining the Si . The main piece
of information which has not yet been used is continuity of S (x). From (3.17) and (3.19) we have
at the ith knot,
S (x) =
(x xi+1 )2 hi+1
fi+1 fi
(x xi )2
Si+1
Si
(Si+1 Si ) +
.
2hi+1
2hi+1
6
hi+1
(3.24)
(x xi1 )2 (x xi )2
hi
fi fi1
Si
Si1 (Si Si1
)+
.
2hi
2hi
6
hi
(3.25)
At x = xi , these two expressions must have the same value, by continuity. Thus, after some
manipulation, we obtain
1
hi+1
fi+1 fi fi fi1
hi
S
+ (hi+1 + hi )Si +
Si+1 =
,
6 i1 3
6
hi+1
hi
(3.26)
which holds for i = 2, . . . , n 1. (It does not hold at the endpoints since only one formula for S is
available at these points.) When i = 1 and x = x1 , we have from (3.24)
(x1 x2 )2 h2
f2 f1
S1
(S2 S1 ) +
2h2
6
h2
h2 h2 f2 f1
S
S +
.
=
3 1
6 2
h2
S (x1 ) =
(3.27)
Sn
(Sn Sn1
)+
2hn
6
hn
hn fn fn1
hn
S
+
S +
.
=
6 n1
3 n
hn
S (xn ) =
(3.28)
58
If (3.27) and (3.28) are used to complete the system (3.26), then endpoint values of S must be
prescribed. In some cases, they may actually be given as part of the data, but it is more usual to
approximate them from the fi s near the endpoints. Recall from our error analysis that we want to
obtain S as a second-order approximation to f in order to maintain fourth-order accuracy for f .
We also found that S must be at least third-order accurate. In a later section we will treat such
derivative approximations. For the present, we will merely denote these required derivatives by f1
and fn , respectively, and assume they are of the required accuracy. Then the system of equations
to be solved for the moments of S takes the form
6 f2 f1
f1 , for i = 1 ;
2S1 + S2 =
h2
h2
hi
fi+1 fi fi fi1
hi+1
6
S
+ 2Si +
S
=
hi + hi+1 i1
hi + hi+1 i+1 hi + hi+1
hi+1
hi
for i = 2, . . . , n 1; and
Sn1
+ 2Sn =
If we define
i
and
6
hn
fn
hi
,
hi + hi+1
fn fn1
hn
i
hi+1
,
hi + hi+1
6 f2 f1
f1
h2
h2
6
fi+1 fi fi fi1
bi =
hi + hi+1
hi+1
hi
fn fn1
6
n
hn
hn
2 1
2 2
for i = n.
for i = 1 ,
for i = n ,
the form
2
2
..
.
3
..
.
n1
S1
b1
S b
0
2 2
= .
..
.
2 n1
S
b
n
n
n
2
(3.29)
This is a tridiagonal system that can be readily solved using the tridiagonal LU decomposition
discussed in Chapter 1.
We point out that there are two other commonly used sets of conditions for determining the
moments of S. The first is simply S1 = Sn = 0. This condition yields what is often called
the natural cubic spline. Observe that this implies that there must be no curvature in f (x) near
the endpoints of the interval [a, b]. It is a widely used alternative, but such wide application is
due mainly to its convenience. It cannot be generally recommended. When it is employed the
system (3.29) contains two fewer equations; namely equations (3.27) and (3.28) are deleted. The
tridiagonal character of the system is not altered because for i = 2, the first term in (3.26) is zero
59
and this is also true for the last term when i = n 1. The final condition sometimes employed to
augment (3.26) is a periodicity condition. It is not widely used, and we will not discuss it further
here. Details can be found in DeBoor [5], a basic practical reference for splines, and in Stoer and
Bulirsch [30].
We shall close this section with a complete algorithm for construction and evaluation of cubic
splines.
Algorithm 3.2 (Cubic Spline with 1st Derivative Specified at Endpoints)
A. Spline Construction
1. Do i = 1, n 1
hi+1 = xi+1 xi
2. Do i = 1, n
a2,i = 2
If i = 1, then a1,i = 0 ;
a3,i = 1 ;
f2 f1
f1
h2
hi
If 1 < i < n, then a1,i = i =
hi + hi+1
hi+1
a3,i = i =
hi + hi+1
6
fi+1 fi fi fi1
bi =
hi + hi+1
hi+1
hi
If i = n then a1,i = 1
6
bi =
h2
a3,i = 0
6
fn fn1
bi =
fn
hn
hn
3. Call tridiagonal LU decomposition to calculate Si , i = 1, . . . , n.
4. Calculate local cubic polynomial coefficients
Do i = 1, n 1
C1,i = fi
fi+1 fi 1
(2Si + Si+1
)hi+1
C2,i =
hi+1
6
1
C3,i = Si
2
1 Si+1 Si
C4,i =
6 hi+1
60
3. Return
3.1.4
f = f + Cj,i (x xi )j1
Extraplotation
3.2
Numerical Quadrature
The term numerical quadrature refers to the numerical evaluation of definite integrals, for example,
of the form
Z
b
f (x) dx .
In general, either or both of a and b may be transfinite. However in our treatment here we will be
concerned only with intervals [a, b] that are of finite length.
The underlying idea in all numerical quadrature is to replace the integrand, f (x), for which
we presumably cannot find a primitive (antiderivative), with an approximation that can be easily
integrated. One of the main approaches to this goal is to approximate f with Lagrange polynomials,
pn , developed in the preceding section, and then integrate pn . This leads to the NewtonCotes
quadrature formulas, of which trapezoidal and Simpsons rules are familiar examples. Another
possibility is to approximate f with Legendre polynomials (which we have not discussed, see, e.g.
[30]), and obtain the GaussLegendre quadrature methods. In any case quadrature formulas can
always be written in the form
Z b
nw
X
wi fi ,
f (x) dx
=h
a
i=1
where the fi are given values of f (x) as appear, for example, in Lagrange interpolation; h is the step
size (usuallybut not alwaysthe distance between successive abscissas where f is evaluated), and
the wi are the nw quadrature weights. It is important to note that the wi vary from one method to
another, but that they do not depend on the function f being integrated. This makes numerical
evaluation of definite integrals extremely simple and efficient, and even very accurate methods can
be implemented on hand-held calculators.
We will here consider three different methods: i) trapezoidal rule, ii) Simpsons rule, and iii)
GaussLegendre quadrature. In addition, we will discuss some methods for improving the accuracy
of trapezoidal quadrature, and for applying any given method to evaluation of multiple integrals.
3.2.1
In this subsection we will consider the two most fundamental NewtonCotes formulas, the trapezoidal rule, and Simpsons rule. We will also discuss two modifications of the former that significantly improve its accuracy.
61
Error at
ith panel
f ( x)
f
fi
h
a
xi+1
xi
n1
b
a
f (x) dx
hX
(fi + fi+1 ) .
2
i=1
Now observe that, except for f1 and fn , each value of f appears twice. Hence, the above can
be written as
#
"
Z b
n1
X
1
fi .
f (x) dx h (f1 + fn ) +
2
a
i=2
We see for this composite (covering more than one panel) trapezoidal formula that the weights are
given by
(
1
if i = 1 or i = n ,
wi = 2
1 otherwise .
It is of interest to note that these weights have the favorable properties, with respect to rounding
errors, of being all of the same sign and being of approximately the same magnitude. Many
quadrature formulas do not possess these attributes.
62
We now present an heuristic argument to show that the trapezoidal rule is second-order accurate;
that is, the error of this approximation is O(h2 ). Recall that we observed at the beginning that
quadrature schemes are generally constructed by integrating polynomial approximations to the
integrand. The above development depends on this, but has been carried out in a very intuitive
manner. Here we will make this slightly more formal, but not completely rigorous. For a complete,
rigorous treatment the reader is referred to the standard text on numerical quadrature, Davis and
Rabinowitz [4].
From Fig. 3.5 it is clear that we have replaced f (x) with a piecewise linear function; so for
x [xi , xi+1 ] [a, b], we have
pi1 (x) = f (x) + O (x xi )2 ,
from results in Section 3.1.2. Now integrate f over the ith panel to obtain
Z xi+1
Z xi+1
pi1 (x)dx + O (xi+1 xi )3 ,
f (x) dx =
xi
xi
and note that O (xi+1 xi )3 O(h3 ). By a basic theorem from integral calculus we have for the
integral over the entire interval [a, b],
!
Z b
n1
n1
n1
X
X Z xi+1
X Z xi+1
h3 .
pi1 (x) dx + O
f (x) dx =
f (x) dx =
a
i=1
xi
i=1
xi
i=1
The first term on the far right is simply the trapezoidal approximation, while for the second term
we have
n1
X
h3 = (n 1)h3 .
i=1
It follows that
Z
b
a
h3 = (b a)h2 .
#
n1
X
1
fi + O(h2 ) .
f (x) dx = h (f1 + fn ) +
2
"
(3.30)
i=2
For the interested reader, we note that it can be shown (cf. [4], or Henrici [13]) that the dominant
truncation error for the trapezoidal method, here denoted simply by O(h2 ), is actually of the form
h2
63
For trapezoidal quadrature it turns out that the exact truncation error on any given subinterval
[xi , xi+1 ] [a, b] can be found from a completely alternative derivation. If instead of employing
Lagrange polynomials to approximate f we use Hermite polynomials (which, again, we have not discussed),
we obtain a fourth-order approximation of f , and thus a locally fifth-order approximation
R
to f . (For details see [30].) In particular, we have
Z xi+1
h3
h5 (4)
h
fi )
f () + ,
(3.31)
f (x) dx = (fi + fi+1 ) (fi+1
2
12
720
xi
provided f C 4 (xi , xi+1 ). Now observe that the first term on the right is exactly the original local
4
fi (fn f1 ) + O(h4 ) .
(3.32)
f (x) dx = h (f1 + fn ) +
2
12
a
i=2
This is called the trapezoidal rule with end correction because the additional terms contain information only from the ends of the interval of integration. If f1 and fn are known exactly, or can be
approximated to at least third-order in h, then (3.32) is a fourth-order accurate method.
The second modification to the trapezoidal rule involves use of extrapolation to cancel leading
terms in the truncation error expansion. A procedure by means of which this can be accomplished
is known as Richardson extrapolation, and this can be used in any situation in which i) the domain
on which the approximations are being done is discretized, and ii) an asymptotic expansion of the
error in powers of the discretization step size is known. In contrast to endpoint correction, the error
terms need not be known exactly for application of Richardson extrapolation. Only the power of
the discretization step size in each term is required.
For the trapezoidal rule, it can be shown (see [30]) that
T (h) = I + 1 h2 + 2 h4 + + m h2m + ,
where I is the exact value of the integral,
I=
and
f (x) dx ,
a
#
n1
X
1
fi ,
T (h) = h (f1 + fn ) +
2
"
(3.33)
i=2
64
Now observe that the dominant error term in (3.34) is exactly 1/4 that in (3.33) since both expansions contain the same coefficients, i . Thus, without having to know the i we can eliminate this
dominant error by multiplying (3.34) by four, and substracting (3.33) to obtain
3
h
T (h) = 3I 2 h4 + O(h6 ) .
4T
2
4
Then division by three leads to the new estimate of I which is accurate to fourth-order:
T (h)
4T ( h2 ) T (h)
1
= I 2 h4 + O(h6 ) .
3
4
(3.35)
An important point to note here is that not only has the original dominant truncation error been
removed completely, but in addition the new dominant term has a coefficient only 1/4 the size of
the corresponding term in the original expansion.
When this procedure is applied recursively to the trapezoidal rule, two orders of accuracy
are gained with each application. This occurs because only even powers of h occur in the error
expansion, as can be seen from (3.33). This technique can be implemented as an automatic highlyefficient procedure for approximating definite integrals known as Romberg integration. Details are
given in [30], and elsewhere.
Simpsons Rule Quadrature
We now briefly treat Simpsons rule. There are several ways to derive this fourth-order quadrature
method. The basic theoretical approach is to replace the integrand of the required integral with a
Lagrange cubic polynomial, and integrate. In Hornbeck [14] Simpsons rule is obtained by a Taylor
expansion of an associated indefinite integral. Here, we will use Richardson extrapolation applied
to the trapezoidal formula. To do this we must employ the global form of the trapezoidal rule,
Eq. (3.30), because we wish to exploit a useful property of the truncation error expansion. As we
have already seen, the global trapezoidal rule has an even-power error expansion while the local
formula contains all (integer) powers of h greater than the second. Thus, recalling Eq. (3.30), we
have
#
"
Z b
n1
X
1
(3.36)
fi + O(h2 ) .
f (x) dx = h (f1 + fn ) +
2
a
i=2
Now we double the value of h and observe (see Fig. 3.6) that on the resulting new partition of
[a, b] only the odd-indexed points of the original partition still occur. In particular, the summation
(3.36) must now run over only odd integers from i = 3 to n 2. This implies that n 2 must be
odd, and hence n is odd. We then have
Z b
n2
X
1
f dx = 2h (f1 + fn ) +
(3.37)
fi + O(h2 ) .
2
a
i=3
i, odd
We now apply Richardson extrapolation by multiplying Eq. (3.36) by four, substracting Eq. (3.37)
and dividing by three:
#
"
Z b
n1
n2
X
X
h
1
1
fi 2 (f1 + fn ) +
f dx =
fi + O(h4 ) .
4 (f1 + fn ) +
3
2
2
a
i=2
i=3
i, odd
65
a
1
.........
i 1 i i+1
h
.......
n 2 n1 n
2h
a
1
b
3
.......
..........
n 2
We can rearrange this expression to obtain a more convenient form by observing that the first sum
contains both even- and odd-indexed terms. Thus,
Z b
n2
n1
X
X
h
fi + 4
fi + O(h4 ) .
f dx = f1 + fn + 2
3
a
i=3
i, odd
i=2
i, even
Finally, it is common for purposes of computer implementation to re-index and write this as
n1
n1
Z b
2
2
X
X
h
f2i1 + 4
f2i + O(h4 ) .
(3.38)
f (x) dx =
f1 + fn + 2
3
a
i=2
i=1
This is the form of Simpsons rule from which very efficient algorithms can be constructed. We see
from (3.38) that the weights for composite Simpsons rule are as follows:
3 for i = 1 or i = n
wi = 23 for 3 i n 2, i odd
4
for 2 i n 1, i even .
3
Also observe that (3.38), as well as the formula that precedes it, reduces to the familiar local
Simpsons rule when n = 3.
3.2.2
GaussLegendre quadrature
The two NewtonCotes methods considered in the preceding section require that function values
be known at equally spaced points, including the endpoints of integration. It turns out that higherorder methods can be constructed using the same number of function evaluations if the abscissas
are not equally spaced. The GaussLegendre quadrature scheme to be considered now is a case of
this. In particular, a GaussLegendre formula employing only n abscissas has essentially the same
accuracy as a NewtonCotes formula using 2n 1 points. Thus, only about half as many integrand
evaluations are required by GaussLegendre to achieve accuracy equivalent to a NewtonCotes
quadrature method. However, this sort of comparison is not very precise and should be viewed as
66
providing only a rough guideline. The more precise statement is the following: A GaussLegendre
method using n absissas will exactly integrate a polynomial of degree 2n 1. By contrast a local
NewtonCotes formula requiring n points will exactly integrate a polynomial of degree n 1.
The GaussLegendre formulas are always local in the sense that the interval of integration is
always [1, 1]. This is because the Legendre polynomials from which the methods are derived
are defined only on [1, 1]. This, however, is not really a serious limitation because any interval
[a, b] R1 can be mapped to [1, 1]. We will here consider only the case of finite [a, b]. If we map
a to 1 and b to 1 by a linear mapping we have
y (1)
1 (1)
2
=
=
,
xa
ba
ba
where x [a, b] and y [1, 1]. Thus, given y [1, 1], we can find x [a, b] from
1
x = a + (b a)(y + 1) .
2
Now, suppose we wish to evaluate
Z
(3.39)
f (x) dx
1
f (x) dx = (b a)
2
f
1
1
a + (b a)(y + 1) dy ,
2
(3.40)
f (x) dx = h
n
X
wi fi .
i=1
For the NewtonCotes formulas h was always the uniform partition step size; but for Gauss
Legendre there is no corresponding quantity. However, if we recall the form of the transformed
interval given above, we see that
ba
.
h=
2
As can be inferred from our discussion to this point the fi do not correspond to evaluations of f
at points of a uniform partition of [1, 1]. Instead the fi are obtained as f (yi ) where the yi are
the zeros of the Legendre polynomial of degree n + 1. Tables of the yi and wi are to be found, for
example in Davis and Rabinowitz [4]. We provide an abbreviated table below for n = 1, 2 and 3.
3.2.3
We will conclude our treatment of basic quadrature methods with a brief discussion of numerical
evaluation of multiple integrals. A standard reference is Stroud [33]. Any of the methods discussed
above can be easily applied in this case; it should be emphasized, however, that the large number of
67
yi
wi
1
2
0.5773502692
0.0000000000
0.7745966692
0.3399810436
0.8611363116
1.0000000000
0.8888888889
0.5555555556
0.6521451547
0.3478548451
function evaluations generally required of the NewtonCotes formulas often makes them unsuitable
when high accuracy is required for triple (or highrer-order) integrals, although such difficulties can
now be mitigated via parallel processing. Here, we will treat only the case of double integrals, but
the procedure employed is easily extended to integrals over domains of dimension higher than two.
Consider evaluation of the integral of f (x, y) over the Cartesian product [a, b] [c, d]. It is
not necessary to restrict our methods to the rectangular case, but this is simplest for purposes of
demonstration. Moreover, nonrectangular domains can always be transformed to rectangular ones
by a suitable change of coordinates. Thus, we evaluate
Z bZ d
f (x, y) dydx .
a
If we define
g(x)
f (x, y) dy ,
(3.41)
we see that evaluation of the double integral reduces to evaluation of a sequence of single integrals.
In particular, we have
Z b
Z bZ d
g(x) dx ;
f (x, y) dydx =
a
so if we set
g(x) dx = hx
a
wi gi ,
(3.42)
i=1
m
X
f (xi , y) dy = hy
c
n
X
wj fij .
j=1
(3.43)
i,j=1
All that is necessary is to choose partitions of [a, b] and [c, d] to obtain hx and hy (unless Gauss
Legendre quadrature is used for one, or both, intervals), and then select a methodwhich determines the wi and wj . We note that it is not necessary to use the same method in each direction,
although this is typically done. We also note that in the context of implementations on modern
parallel processors, it is far more efficient to evaluate the m equations of the form (3.41) in parallel,
and then evaluate (3.42) instead of using (3.43) directly.
68
3.3
Finite-Difference Approximations
Approximation of derivatives is one of the most important and widely-used techniques in numerical
analysis, mainly because numerical methods represent the only general approach to the solution of
differential equationsthe topic to be treated in the final two chapters of these lectures. In this
section we will present a formal discussion of difference approximations to differential operators.
We begin with a basic approximation obtained from the definition of the derivative. We then
demonstrate use of Taylor series to derive derivative approximations, and analyze their accuracy.
Following this we will consider approximation of partial derivatives and derivatives of higher order.
We then conclude the section with a few remarks and approximation methods that are somewhat
different from, but still related to, the finite-difference approximations described here.
3.3.1
Basic concepts
We have already used some straightforward difference approximations in constructing the secant
method and cubic spline interpolation. These basic approximations follow from the definition of
the derivative, as given in Freshman calculus:
lim
h0
f (x + h) f (x)
= f (x) ,
h
(3.44)
provided the limit exists. To obtain a finite-difference approximation to f (x) we simply delete the
limit operation. The result is the first forward difference,
f (x + h) f (x)
.
h
= xi + h, then in our usual notation we see that
f (x)
If we note that on a grid of points xi+1
fi+1 fi
fi+1 fi
=
h
xi+1 xi
(3.45)
1
fi + fi h + fi h2 +
2
1
= fi + fi h + .
2
fi
Hence, the leading error in (3.45) is 21 fi h; so the approximation is first order in the step size h.
3.3.2
There are many different ways to obtain derivative approximations, but probably the most common is by means of the Taylor series. We will demonstrate this now for a backward-difference
approximation to the first derivative. We again assume f C 2 , and write
1
fi1 = fi fi h + fi h2 .
2
69
fi fi1
+ O(h) .
h
(3.46)
In order to obtain derivative approximations of higher-order accuracy, we can carry out Taylor
expansions to higher order and form linear combinations so as to eliminate error terms at the
desired order(s). For example, we have
1
1
1
1 5
fi+1 = fi + fi h + fi h2 + fi h3 + fi h4 +
f h + ,
2
6
24
120 i
and
1
1
1
1 5
fi1 = fi fi h + fi h2 fi h3 + fi h4
f h .
2
6
24
120 i
If we substract the second from the first, we obtain
1
fi+1 fi1 = 2fi h + fi h3 + ,
3
and division by 2h leads to
fi =
fi+1 fi1 1 2
1 4
fi h
f h .
2h
6
120 i
(3.47)
+ O(h4 )
fi =
3
2h
4h
(3.48)
1
4
(fi2 8fi1 8fi+1 fi+2 ) + O(h ) .
=
12h
This is the fourth-order accurate centered approximation to the first derivative.
There is yet another way to employ a Taylor expansion of a function to obtain a higher-order
difference approximation. This involves expressing derivatives in the leading truncation error terms
as low-order difference approximations to eliminate them in favor of (known) grid function values.
We demonstrate this by constructing a second-order forward approximation to the first derivative.
Since we are deriving a forward approximation we expect to use the value of f at xi+1 . Thus, we
begin with
1
1
fi+1 = fi + fi h + fi h2 + fi h3 + ,
2
6
and rearrange this as
fi+1 fi 12 fi h2 1 2
fi =
fi h + .
(3.49)
h
6
We now observe that we can obtain an approximation of the desired order if we have merely a
first-order approximation to fi . We have not yet discussed approximation of higher derivatives,
70
but as we will see later the only required idea is simply to mimic what we do analytically for exact
derivatives; namely we repeatedly apply the difference approximation. Now recall that
fi =
fi+1 fi
+ O(h) ;
h
fi =
fi+1
fi
+ O(h) .
h
fi =
fi =
fi+1 fi (fi+2 2fi+1 + fi ) + O(h2 ) ,
h
2
or, after rearrangement,
1
(3fi + 4fi+1 fi+2 ) + O(h2 ) .
(3.50)
2h
This is the desired second-order forward approximation. A completely analogous treatment leads
to the second-order backward approximation; this is left as an exercise for the reader.
fi =
3.3.3
The next topics to be treated in this section are approximation of partial derivatives and approximation of higher-order derivatives. We will use this opportunity to introduce some formal notation
that is particularly helpful in developing approximations to high-order derivatives and, as will be
seen in the next two chapters, for providing concise formulas for discrete approximations to differential equations. The notation for difference approximations varies widely, and that used here is
simply the preference of this author.
In general we will take D(h) to be a difference operator based on step size h. (When no confusion
is possible we will suppress the notation for h.) We then denote the forward-difference operator
by D+ (h), the backward operator by D (h) and centered operators by D0 (h). Thus, we have the
following:
fi+1 fi
= f (xi ) + O(h) ,
h
fi fi1
D (h)fi =
= f (xi ) + O(h) ,
h
fi+1 fi1
= f (xi ) + O(h2 ) .
D0 (h)fi =
2h
D+ (h)fi =
(forward)
(3.51a)
(backward)
(3.51b)
(centered)
(3.51c)
When we require partial derivative approximations, say of a function f (x, y), we alter the above
notation appropriately with, for example, either Dx (h) or Dy (h). Hence, for the centered difference
we have
fi+1,j fi1,j
f
=
(xi , yi ) + O(h2 ) ,
(3.52)
D0,x (h)fi,j =
2h
x
f
fi,j+1 fi,j1
=
(xi , yi ) + O(h2 ) .
2h
y
71
(3.53)
We noted earlier that approximation of higher derivatives is carried out in a manner completely
analogous to what is done in deriving analytical derivative formulas. Namely, we utilize the fact
that the (n + 1)th derivative is just the (first) derivative of the nth derivative:
dn+1
d dn f
f=
.
dxn+1
dx dxn
In particular, in difference-operator notation, we have
Dn+1 (h)fi = D(h) (D n (h)fi ) .
We previously used this to obtain a first-order approximation of f , but without the formal notation.
We will now derive the centered second-order approximation:
1
2
D0 fi = D0 (D0 fi ) = D0
(fi+1 fi1 )
2h
1
fi+2 fi
fi fi2
=
2h
2h
2h
1
(fi2 2fi + fi+2 ) .
=
(2h)2
We observe that this approximation exactly corresponds to a step size of 2h, rather than to h,
since all indices are incremented by two, and only 2h appears explicitly. Hence, it is clear that the
approximation over a step size h is
D02 fi =
(3.54)
In recursive construction of centered schemes, approximations containing more than the required
range of grid point indices always occur because the basic centered operator spans a distance 2h. It
is left to the reader to verify that (3.54) can be obtained directly, by using the appropriate definition
of D0 (h) in terms of indices i 21 and i + 12 . We also note that it is more common to derive this
using a combination of forward and backward first-order differences, D+ D fi .
3.3.4
There are two remaining approximation methods which are related to differencing, and which are
widely used. The first is divided differences. We will not treat this method here, but instead discuss
the second approach, which gives identical results. It is simply differentiation of the Lagrange (or
other) interpolation polynomial. Suppose we are required to produce a second-order accurate
derivative approximation at the point xi . Now we expect, on the basis of earlier discussions, that
differentiation of a polynomial approximation will reduce the order of accuracy by one power of the
step size. Thus, if we need a first derivative approximation that is second-order accurate, we must
start with a polynomial which approximates functions to third order.
Hence, we require a quadratic, which we formally express as
p2 (x) =
3
X
i=1
72
3
X
i=1
i (x)fi + O(h2 ) .
We can now obtain values of f for any x [x1 , x3 ]. In general, we typically choose the x so that
x = xi for some i. The main advantage of this Lagrange polynomial approach is that it does not
require uniform spacing of the xi s, such as is required by all of the procedures presented earlier. (It
should be noted, however, that Taylor series methods can also be employed to develop difference
approximations over nonequally spaced points; but we shall not pursue this here).
3.4
We have previously used Richardson extrapolation to construct Simpsons rule from trapezoidal
quadrature, and also to obtain higher-order difference approximations; it is also the basis for the
Romberg integration method. In recent years Richardson extrapolation has come into wide use in
several important areas of computational numerical analysis, and because of this we feel a more
general treatment is needed than can be deduced merely from the specific applications discussed
above. In all of these examples we were extrapolating a procedure from second- to fourth-order
accuracy, which depends upon the fact that only even powers of the step size appear in the truncation error expansion. Furthermore, extrapolation was always done between step sizes differing by a
factor of two. There are many cases in which all powers of h (and for that matter, not just integer
powers) may occur in the truncation error expansion. Moreover, for one reason or another, it may
not be convenient (or even possible) to employ step sizes differing by a factor of two. Hence, it is
important to be able to construct the extrapolation procedure in a general way so as to remove
these two restrictions.
Let {xi }ni=1 be a partition of the interval [a, b] corresponding to a uniform step size h = xi+1 xi .
Let f (xi ) be the exact values of a function f defined on this interval, and let {fih }ni=1 denote the
corresponding numerical approximation. Hence,
fih = f (xi ) + 1 hq1 + 2 hq2 + O (hq3 )
(3.55)
for some known qm R, m = 1, 2, . . . , and (possibly) unknown m C which also depend on the
grid point xi . We have earlier seen in a special case that it is not necessary to know the m .
The functions fih , usually called grid functions, may arise in essentially any way. They may
result from interpolation or differencing of f ; they may be a definite integral of some other function,
say g, as in our quadrature formulas discussed earlier (in which case the i index is superfluous); or
they might be the approximate solution to some differential or integral equation. We will not here
need to be concerned with the origin of the fih s.
Let us now suppose that a second approximation to f (x) has been obtained on a partition of
[a, b] with spacing rh, r > 0 and r 6= 1. We represent this as
firh = f (xi ) + 1 (rh)q1 + 2 (rh)q2 + O (hq3 ) .
(3.56)
We note here that we must suppose there are points xi common to the two partitions of [a, b].
Clearly, this is not a serious restriction when r is an integer or the reciprocal of an integer. In
general, if the fi s are being produced by any operation other than interpolation, we can always, in
73
principle, employ high-order interpolation formulas to guarantee that the grid function values are
known at common values of x [a, b] for both values of step size.
We now rewrite (3.56) as
firh = f (xi ) + r q1 1 hq1 + r q2 2 hq2 + O (hq3 ) .
From this it is clear that the q1th -order error term can be removed by multiplying (3.55) by r q1 , and
substracting this from (3.56). This leads to
firh r q1 fih = f (xi ) r q1 f (xi ) + r q2 2 hq2 r q1 2 hq2 + O (hq3 )
= (1 r qi ) f (xi ) + O (hq2 ) .
From this it is clear that we should define the general extrapolated quantity fi as
fi
r q1 fih firh
= f (xi ) + O (hq2 ) .
r q1 1
(3.57)
We now demonstrate for the special cases treated earlier, for which q1 = 2, q2 = 4, and r = 2,
that the same result is obtained using Eq. (3.57) as found previously; namely
fi =
1 h
4fih f 2h
4fi fi2h = f (xi ) + O(h4 ) .
=
41
3
Similarly, if the leading term in the truncation error expansion is first order, as was the case with
several of the difference approximations presented in the preceding section, we have q1 = 1, q2 = 2,
and for r = 2 (3.57) yields
fi = 2fih fi2h = f (xi ) + O(h2 ) .
We again emphasize that the Richardson extrapolation procedure can be applied to either
computed numerical results (numbers), or to the discrete formulas used to produce the results.
The optimal choice between these alternatives is typically somewhat problem dependentthere
are no general prescriptions.
3.5
74
the theoretical convergence rate of the grid functions generated by the method. In particular, we
almost always have the truncation error expansion at our disposal. We have derived numerous of
these throughout this chapter.
For example, from (3.55) we see that
fih = f (xi ) + 1 hq1 + = f (xi ) + O (hq1 ) ,
and by changing the step size to rh we have
firh = f (xi ) + 1 r q1 hq1 + .
The dominant error in the first case is
ehi f (xi ) fih = 1 hq1 ,
(3.58)
rh
q1 q1
erh
,
i = f (xi ) fi = 1 r h
(3.59)
(3.60)
Hence, for a second-order method (q1 = 2) a reduction in the step size by a factor of two (r = 12 )
leads to a reduction in error given by
2
1
1
q1
= ;
r =
2
4
i.e., the error is reduced by a factor of four.
In practical problems we usually do not know the exact solution, f (x); hencenwe cannot
calculate
o
o
n
h
h/4
h/2
,
and fi
the true error. However, if we obtain three approximations to f (x), say fi , fi
we can make good estimates of 1 , q1 and f (xi ) at all points xi for which elements of all three grid
functions are available. This merely involves solving the following system of three equations for 1 ,
q1 and f (xi ):
fih = f (xi ) + 1 hq1 ,
h/2
h/4
fi
fi
h/2
h/4
h/2
fi
h/4
fi
= 2q1 1 2q1 1 hq1 .
(3.62)
= 2q1 ,
(3.63)
h/2
fih fi
h/2
fi
h/4
fi
75
which is equivalent to the result (3.60) obtained above using true error. Again note that q1 should
be known, theoretically; but in practice, due either to algorithm/coding errors or simply to use of
step sizes that are too large, the theoretical value of q1 may not be attained at all (or possibly at
any!) grid points xi .
This motivates us to solve Eq. (3.63) for the actual value of q1 :
#
"
h/2
fih fi
log h/2
h/4
fi fi
.
(3.64)
q1 =
log 2
Then from Eq. (3.61) we obtain
h/2
fih fi
1 =
.
(3.65)
1 2qi hq1
Finally, we can now produce an even more accurate estimate of the exact solution (equivalent to
Richardson extrapolation) from any of the original equations; e.g.,
f (xi ) = fih 1 hq1 .
(3.66)
In most practical situations we are more interested in simply determining whether the grid
functions converge and, if so, whether convergence is at the expected theoretical rate. To do this it
is usually sufficient to replace f (xi ) in the original expansions with a value fi computed on a grid
much finer than any of the test grids, or a Richardson extrapolated value obtained from the test
grids, say fi . The latter is clearly more practical, and for sufficiently small h it leads to
1 hq1 .
eh = f f h =
i
Similarly,
h/2
ei
and the ratio of these errors is
h/2
= fi fi
= 2q1 1 hq1 ,
ehi fi fih
= 2q1 .
(3.67)
=
h/2
h/2
ei
fi fi
Yet another alternative (and in general, probably the best one when only grid function convergence is the concern) is simply to use Eq. (3.63), i.e., employ a Cauchy convergence test. As noted
above we generally know the theoretical value of q1 . Thus, the left side (obtained from numerical
computation) can be compared with the right side (theoretical). Even when q1 is not known we
can gain qualitative information from the left-hand side alone. In particular, it is clear that the
right-hand side is always greater than unity. Hence, this should be true of the left-hand side. If the
equality in the appropriate one of (3.63) or (3.67) is not at least approximately satisfied, the first
thing to do is reduce h, and repeat the analysis. If this does not lead to closer agreement between
left- and right-hand sides in these formulas, it is fairly certain that there are errors in the algorithm
and/or its implementation.
We note that the above procedures can be carried out for arbitrary sequences of grid spacings,
and for multi-dimensional grid functions. But in both cases the required formulas are more involved,
and we leave investigation of these ideas as exercises for the reader. Finally, we must recognize
that ehi (or ehi ) is error at a single grid point. In most practical problems it is more appropriate to
employ an error norm computed with the entire solution vector. Then (3.63), for example, would
be replaced with
keh k kf h f h/2 k
= 2q1 ,
(3.68)
= h/2
keh/2 k
kf
f h/4 k
for some norm k k, say, the vector 2-norm.
76
3.6
Summary
This chapter has been devoted to presenting a series of topics which, taken together, might be called
classical numerical analysis. They often comprise a first course in numerical analysis consisting
of interpolation, quadrature and divided differences. As will be evident in the sequel, these topics
provide the main tools for development of numerical methods for differential equations.
As we have attempted to do throughout these lectures, we have limited the material of this
chapter to only the most basic methods from each class. But we wish to emphasize that, indeed,
these are also the most widely used for practical problem solving. All of the algorithms presented
herein (and many similar ones) are available in various commercial software suites, and above all
else we hope the discussions presented here will have provided the reader with some intuition into
the workings of such software.
Chapter 4
4.1
Initial-Value Problems
In this section we discuss some of the main methods for solving initial-value problems for ordinary
differential equations. We begin with a brief mathematical background for such problems, and
then proceed to treat single-step numerical methods, multi-step methods and finally techniques for
solving stiff ODEs.
4.1.1
Mathematical Background
Initial-value problems for nth -order ODEs can always be cast in the following general form:
F u(n) , u(n1) , . . . , u , u, t = 0
(4.1)
77
78
(4.2)
u(n1) (t0 ) = cn ,
where F, u Rm and t R1 . Parenthesized superscripts here denote order of differentiation.
We should start by noting that an nth -order differential equation will require n initial conditions,
as given in (4.2), to permit evaluation of the n integration constants arising during the n formal
integrations needed to obtain a solution.
We will not usually encounter systems of ODEs in a form as general as (4.1), (4.2). It is often,
but not always, possible to solve each equation in the system for its highest derivative term. When
this is the case (4.1) can be written as
dn u
(n1) (n2)
=
f
u
,
u
,
.
.
.
,
u
,
u,
t
,
dtn
(4.3)
with f Rm related to F in an obvious way. The initial data (4.2) still apply in this case.
It is an important fact that any equation of the form (4.3) can be written as a system of n firstorder equations. Moreover, initial conditions of the form (4.2) provide the correct initial data for
this system, as we will show below. The consequence of all this is that study of the numerical initialvalue problem can be restricted to first-order systems without loss of generality. This approach,
although not the only possibility, is essentially always employed, and we will follow it here.
To express (4.3) as a first-order system, let
du
y1 =
,
dt
d2 u
= 2
dt
dy1
y2 =
dt
yn1
dyn2
=
dt
dn1 u
= n1
dt
dyn2
dt
d3 yn3
dn1 y1
dn u
d2 yn2
=
=
=
=
=f ,
dt2
dt3
dtn1
dtn
dyn2
= yn1
dt
dyn1
= f (yn1 , yn2 , . . . , y1 , u, t) .
dt
(4.4)
79
This is a system of n equations for u, the desired solution, and its first n 1 derivatives. The
required initial data are given in (4.2). In particular, we have
u(t0 ) = c1
u (t0 ) = y1 (t0 ) = c2
(4.5)
2
1
du1
du2 2
d2 u1
+ sin
u2 cos t = 0 ,
dt2
dt
dt
2
d u2
du1
d3 u2
exp
+ u1 = 0 .
+ u2
dt3
dt2
dt
(4.6a)
(4.6b)
Since (4.6a) is second order, two initial conditions will be needed, while (4.6b) is third order, and
will require three conditions. In particular, we have
u1 (t0 ) = a
u2 (t0 ) = c
u1 (t0 )
u2 (t0 ) = d
=b
u2 (t0 )
(4.7)
=e.
The first step is to solve each of (4.6a) and (4.6b) for their respective highest derivative. For
(4.6a) we obtain
d2 u1
=
dt2
cos t + u2 +
du2
dt
1
2
sin
du1
dt
)21
f1 ,
(4.8a)
d2 u2
dt2
u2 +
du1
u1 f 2 .
dt
(4.8b)
du1
,
dt
dy1
d2 u1
=
= f1 ,
dt
dt2
(4.9)
z1 =
du2
,
dt
z2 =
dz1
,
dt
dz2
=
dt
d3 u2
dt3
= f2 .
80
= y1
= {cos t + u2 +
u1 (t0 ) = a
1
z1 sin y1 } 2
y1 (t0 ) = b
= z1
u2 (t0 ) = c
= z2
z1 (t0 ) = d
= exp(z2 ) u2 + y1 u1
z2 (t0 ) = e .
(4.10)
The system (4.10) is coupled and strongly nonlinear. Existence and/or uniqueness of solutions for
such systems is not guaranteed, a priori. We shall not address such questions in this presentation;
but the reader must be warned that since numerical methods for initial-value problems can be
applied irrespective of nonlinearities, computed results to problems such as (4.10) must be viewed
with some caution.
A second remark is also important at this time. It is that although we have demonstrated that
any higher-order system of the form (4.3) can be written as an equivalent first-order system, it is
not true that every first-order system arises as a reduction of some higher-order equation. Thus,
any general algorithm for solving ODE IVPs must be written for general systems, and not merely
for the special form shown in (4.4).
4.1.2
In this section we will begin with a detailed treatment of the simplest single-step method, the
forward Euler method; we will then consider the backward Euler procedure. Following this we
present analyses of two higher-order methods, the implicit trapezoidal scheme and its explicit
counterpart Heuns method, which is an explicit second-order RungeKutta method.
Explicit (Forward) Euler Method
The first method we consider, Eulers method, is one that should almost never be used in practice
in the context of ODEs (but we will later see that it is often used in solution techniques for partial
81
differential equations). However, it is extremely simple and thus provides a useful pedagogical tool.
It will be evident that the method is easily coded even on a hand-held calculator, and we might
at times consider using it in this context. We will treat only a single first-order equation since the
extension to systems is mainly a matter of notation.
Consider the equation
u = f (u, t) ,
(4.11)
with initial condition u(t0 ) = u0 . If we replace u with a first forward difference and evaluate f at
time tn we obtain
un+1 un
= f (un , tn ) ,
h
or
un+1
(4.12)
= un + hf (un , tn ) .
Clearly if un is known, we can evaluate the right-hand side and thus explicitly calculate un+1 .
This is a characteristic feature of explicit methods: the grid function values at a new time step can
be directly evaluated from values at previous time stepswithout numerical linear algebra and/or
iteration. We also observe that the right-hand side of (4.12) involves information from only a single
previous time step. Hence, Eulers method is an explicit single-step method.
We now investigate the truncation error for Eulers method. Intuitively, we would expect this
method to be only first-order accurate because we have used a first-order approximation to u .
Indeed this turns out to be true in a certain sense; but it is important to understand some details.
Our treatment is not intended to be rigorous, but it is correct with respect to the basic notions.
(For a rigorous development, we recommend Gear [10]). In going from the ODE to the difference
approximation, the exact result would be (for sufficiently small h)
un+1 un
+ O(h) = f (un , tn ) .
h
It is clear from this that as h 0, the original differential equation is recovered. Thus, Eulers
method is said to be consistent with the differential equation. Analogous to (4.12) we have
un+1 = un + hf (un , tn ) + O(h2 ) ,
(4.13)
and in this form Eulers method would appear to be second-order accurate. However, the formula
(4.13) advances the solution only a single time step, and during this step an error of order h2 is
incurred. This is called the local truncation error, and it is second order for Eulers method.
Now suppose we wish to solve (4.11) on the interval (0, ], using N time steps. Then the
(uniform) time step is h = /(N 1). From (4.13) we see that after N 1 steps (the number needed
to obtain uN ) we will have accumulated a truncation error equal to (N 1) O(h2 ), analogous to
what we found earlier in analysis of trapezoidal quadrature in Chap. 3. Thus, the global truncation
error for Eulers method is O(h). The order of a method is always taken to be the order of the
global truncation error; hence, Eulers method is first-order accurate just as we originally expected.
We next consider the stability of Eulers method. There are many different definitions employed
in studying stability, but basically, we consider a difference approximation to a given problem to
be stable if the solution to the difference equation does not blow up any faster than the solution
to the differential equation. For the analyses presented herein, we employ the following somewhat
more precise statement of this idea.
Definition 4.1 A method is absolutely stable for a given step size h, and for a given differential
equation, if the change due to a perturbation of size in any mesh value um is no larger than
un , n > m.
82
We note that the perturbations referred to in the definition typically arise from round-off error in
machine computation.
To examine absolute stability of Eulers method we consider a very simple IVP,
u = u ,
u(0) = u0 ,
const. < 0 .
(4.14)
un = (1 + h)un1 = = (1 + h)n u0 .
Observe that this is an exact solution to the difference equation corresponding to Eulers method.
Now suppose u0 is replaced by v0 = u0 + , where is the error in machine representation of
u0 and thus corresponds to a perturbation in the sense of the definition of absolute stability. For
example, if u0 = , || > 0 will hold on any machine, no matter what the word size happens to be
because it will always be finite, and does not have a finite exact representation. After n steps we
have
vn = (1 + h)n (u0 + )
= (1 + h)n u0 + (1 + h)n
= un + (1 + h)n .
Now define the error at time step n to be
zn = vn un = (1 + h)n ,
which represents the growth in time of the perturbation .
Then, taking z0 = , we see that the error satisfies the same difference equation as does the
solution, un , and after n steps the original error will have been amplified to
(1 + h)n .
It follows that in order to guarantee absolute stability for Eulers method we must have
|1 + h| 1 .
(4.15)
83
Unstable
Im h
Unstable
Stable
( 1,0)
Re h
Unstable
Unstable
Eulers method is stable for h [2, 0]. This shows that the method is never absolutely stable
if > 0, and for < 0, but || 1, the step sizes required for stability are extremely small. For
these reasons, Eulers method is seldom used in practical calculations.
It is interesting to observe the effects of instability on the computed solution. For simplicity, we
take u0 = 1 in (4.14), and set = 100. At t = 1, the exact solution is u(1)
= 3.72 1044 0.
From (4.15) we see that to guarantee stability we must choose the step size h so that
|1 100h| 1 ,
or
1 100h 1
and 1 100h 1 .
The first of these implies h 0, which we would expect in any case. From the second we find
h
With h exactly equal to
1
50 ,
1
.
50
Hence, as the time steps proceed u(t) = u0 = 1 which remains bounded, but is completely
wrong.
Next consider h = 0.01. In this case, the amplification factor is zero, and un 0 n > 0. This
is very inaccurate close to t = 0, but asymptotically correct as t . If we set h = 0.001, the
amplification factor is 0.9. To reach t = 1 we need 1000 time steps, so the Eulers method solution
is
u1000 = (0.9)1000 u0
= 1.748 1046 ,
which is at least beginning to resemble the exact solution. Moreover, the computed solution at
t = 0.1 is u100
= 2.66 105 compared with the exact solution u(.1) = 4.5 105 . Hence, the
method produces results that are at least of the correct order of magnitude for short times, and
decreasing h to h = 0.0001 yields reasonably accurate results.
84
h = 0.001
h = 0.02
1
h = 0.01
2
0.0
0.2
0.1
h = 0.1
0.3
0.4
0.5
Finally, we consider h = 0.1; then the amplification factor is 9, and Eulers method becomes
un = (9)n u0 .
To integrate to t = 1, we need 10 steps, and this leads to
u10
= 3.487 109 ,
which is completely ridiculous.
Figure 4.2 summarizes these results. Of particular note is the phenomenon of growing oscillations in the unstable case. Such oscillations provide a tell-tale symptom of instability rather
generally in numerical integration procedures, and their appearance is a signal that the time step
size should be reduced.
We conclude from all this that Eulers method is generally unreliable because of its poor stability
characteristics, and it should be used only with extreme caution.
Implicit (Backward) Euler Method
There is a relatively simple remedy to the stability problems encountered with Eulers method.
Recall that in developing this method we replaced u with a forward-difference approximation. We
now instead use a backward approximation to obtain
un un1
= f (un , tn ) ,
h
or
un = un1 + hf (un , tn ) .
If we translate this forward by one time step the result is analogous to the earlier form of Eulers
method (4.12):
un+1 = un + hf (un+1 , tn+1 ) .
(4.16)
85
This approximation is known as the backward Euler method. It is still a single-step method, but it
is now implicit. In particular, we cannot in general calculate a new time step by merely plugging
in results from a preceding time step. Usually, f will be a nonlinear function of un+1 ; so (4.16) is
a nonlinear equation for un+1 that can be solved by a fixed-point algorithm applied at each time
step. Newtons method is typically used. Thus, we write (4.16) as
F (un+1 ) = un+1 un hf (un+1 , tn+1 ) = 0 ,
(4.17)
un+1
(m)
F un+1
(m)
,
= un+1
(m)
F un+1
(4.18)
f
.
un+1
(4.19)
The truncation error analysis for backward Euler can be carried out analogously to what was done
for forward Euler, and the result is the same; the backward Euler method is first-order accurate.
We leave demonstration of this as an exercise for the reader.
We now consider the stability properties of (4.16) with respect to the IVP (4.14). In this case
f is linear in u, so we have
un+1 = un + hun+1 ,
which can be solved for un+1 without iteration (because of linearity):
un+1 hun+1 = un ,
and
un+1 =
un
.
1 h
(4.20)
Analogous to what we did in analyzing absolute stability of (forward) Eulers method, we find for
backward Euler that
1
u0 .
(4.21)
un =
(1 h)n
Thus, the amplification factor is |(1 h)1 |, and the condition for absolute stability is
1
1,
|1 h|
or
|1 h| 1 .
(4.22)
This is all of the complex h-plane except for the open disc of unit radius centered at (1, 0).
Our first observation is that (4.22) is satisfied 0 (or Re 0 if C), independent
of the step size h. Thus, for the problem considered previously, with = 100, there are no
1
, the amplification factor is 13 ,
restrictions on h for maintaining stability. For example, for h = 50
so
n
1
u0 .
un =
3
86
Hence u50 u(1) = 1.393 1024 0. This is actually much larger than the exact result, but
both are so near zero that in many practical situations the backward Euler result would be quite
1
, we have u100 u(1) = 7.889 1031 which is significantly closer
acceptable. If we choose h = 100
to the exact result. Finally, for h = 0.1, which is very unstable for forward Euler, the backward
1 10
= 3.855 1011 , which is still sufficiently close to zero
Euler result at t = 1 is u10 u(1) = 11
for many practical purposes.
We close this section on Euler single-step methods with the following remarks. We have shown
the Euler methods to be stable, at least for sufficiently small step size h, and we have also indicated
that they are consistent with the differential equation whose solution is being approximated. But
we have not proven that the numerical solutions actually converge to the solution of the differential
equation as h 0. Here, we will simply note that such convergence is guaranteed by the combination of consistency and stability. We will consider this more formally and in greater detail, in the
context of partial differential equations in Chap. 5.
Error
Total Error
Truncation
Error
Step size of
minimum
total error
Round-off
Error
Increasing h
87
The first higher-order method we shall treat here is the trapezoidal rule. Once again, we consider
the first-order differential equation
u = f (u, t)
(4.23)
with the initial condition
u(t0 ) = u0 .
(4.24)
Now the left-hand side is exact, and we approximate the integral on the right with the (local)
trapezoidal rule described in Chap. 3, which results in
un+1 = un +
h
[f (un+1 , tn+1 ) + f (un , tn )] + O(h3 ) .
2
(4.25)
Recall that local truncation error for the trapezoidal rule is O(h3 ), and by an argument analogous
to that used for Eulers method, the global truncation error is O(h2 ). Thus, trapezoidal integration
is a second-order method, and as a consequence we expect truncation error to be reduced by a
factor of four each time the step size h is halved.
An important observation to make regarding (4.25) is that it is implicit, and we must employ an
iteration scheme very similar to that used above for the backward Euler method. As a consequence f
must generally be evaluated for each Newton iteration, making this approach somewhat inefficient.
On the other hand, rather generally, implicit methods have very favorable stability properties, and
this is true for the trapezoidal rule. It can be shown that its region of absolute stability is essentially
all of the left-half complex h-plane, as we will indicate later. However, before we present details
of a trapezoidal integration algorithm, it is of interest to consider an explicit method which can be
easily obtained from the trapezoidal rule.
Heuns Method
Observe that in (4.25) we cannot explicitly evaluate f (un+1 , tn+1 ) because we do not yet know
un+1 , but if we could estimate un+1 with sufficient accuracy using only previous information we
could evaluate (4.25) explicitly. In particular, denote this estimate as un+1 , and rewrite (4.25) as
un+1 = un +
h
f (un+1 , tn+1 ) + f (un , tn ) .
2
(4.26)
The only candidate presently available for explicitly calculating un+1 is Eulers method, so we set
un+1 = un + hf (un , tn ) .
Substitution of this into (4.26) leads to
un+1 = un +
h
[f (un + hf (un , tn ), tn+1 ) + f (un , tn )] .
2
(4.27)
This explicit method is known as Heuns method. We will later see that it is one of the simplest
nontrivial RungeKutta schemes. It can also be viewed as a simple predictor-corrector technique.
These will be discussed in more detail later in the study of multi-step methods.
88
We note here that Heuns method is also globally second order, despite the fact that we have
used a first-order approximation to obtain un+1 . To prove this we need to show that
f (un+1 , tn+1 ) = f (un+1 , tn+1 ) + O(h2 ) .
Then we see from (4.26) that the local error is O(h3 ). To prove the above equality we will assume
f is Lipschitz in u, but we note that this is no more restrictive than what is needed to guarantee
uniqueness of solutions to (4.23) in the first place. If f is Lipschitz we have
f (u , tn+1 ) f (un+1 , tn+1 ) L u un+1
n+1
n+1
= L |un + hf (un , tn ) un+1 |
O(h2 ) .
This holds because the expression in absolute values on the right-hand side is just Eulers method,
which we have already shown in Eq. (4.13) to be locally second order. This provides the desired
result.
Heuns method (4.27) is an easily programmed explicit scheme. It is significantly more accurate
than Eulers method, and it also has somewhat better stability properties when applied to problems
having oscillatory solutions. Nevertheless, it is important to note that there are many other higherorder methods which are more accurate, more stable and fairly efficientbut more difficult to
construct.
Heun/Trapezoidal Algorithm
We will now develop a detailed algorithm that permits solution of general systems of first-order
IVPs by either Heuns method or by the trapezoidal rule. As noted earlier, trapezoidal integration
is implicit, and at each time step the (generally) nonlinear system
Fi (un+1 ) ui,n+1 ui,n
h
[f (ui,n+1 , tn+1 ) + f (ui,n , tn )] = 0
2
i = 1, . . . , p ,
(4.28)
with u, f, F Rp must be solved. This is done via Newtons method, so we must supply the
Jacobian matrix of F as part of the problem input. We show that Ju (F ) can be expressed entirely
in terms of the partial derivatives of f and the step size h. This will permit us to prepare input
based only on the analytical problem being solvedthus not requiring problem-dependent changes
to the discrete equations of the integration procedure.
Let
un+1 = (u1,n+1 , u2,n+1 , . . . , up,n+1 )T
for a system of p first-order equations. Then for the ith component of F , (4.28) leads to
h fi
Fi
= ij
uj,n+1
2 uj,n+1
Thus, we see that
i, j = 1, 2, . . . , p .
(4.29)
h
Ju (f ) ,
2
so in order to input a problem to be solved by the trapezoidal method we need only code the fi s
and the elements of Ju (f ).
In addition, at each time step we need to supply an initial guess for the Newton iterations. Since
we will produce an algorithm that also can employ Heuns method, it is reasonable to use this as
Ju (F ) = I
89
the initial guess for the trapezoidal rule at each time step. This leads to very rapid convergence
of the Newton iterations. But we note that it is not really necessary to use such an elaborately
constructed initial guess. It is usually quite acceptable to employ the result from the previous
time step, or possibly an Eulers method prediction. We now summarize the preceding ideas in a
pseudo-language algorithm.
Algorithm 4.1 (Heun/Trapezoidal Method)
1. Read control flags: neqns, maxitr, nsteps
Read numerical parameters: h, ,
Read initial data: (ui,0 , i = 1, neqns), t0
2. Begin time stepping
Do n = 0, nsteps 1
tn+1 = tn + h
3. Begin Newton iterations for current time step
Do m = 1, maxitr
If m > 1, go to 6
4. Evaluate ui,n+1 for use in Heuns method
Do i = 1, neqns
gi = f (i, un , tn )
ui,n+1 = ui,n + hgi
Repeat i
5. Calculate initial guess for trapezoidal rule from Heuns method
Do i = 1, neqns
(0)
ui,n+1 = ui,n + h2 (gi + f (i, un+1 , tn+1 ))
Repeat i
If maxitr = 1, go to 11 [Heuns Method]
6. Load J(f ) for Newton
iteration
(m1)
Call Jacobn neqns, un+1 , tn+1 , J(f )
(m1)
7. Evaluate F un+1
, J(F )
Do i = 1, neqns
(m1)
(m1)
Fi = ui,n+1 h2 f i, un+1 , tn+1 ui,n + h2 gi
Do j = 1, neqns
J(F )ij = ij h2 J(f )ij
Repeat j
Repeat i
90
(m1)
f = f1 (u, t)
Return
20
f = f2 (u, t)
Return
End
Subroutine Jacobn(neqns, u, t, J(f ))
Do i = 1, neqns
Do j = 1, neqns
J(f )ij = fi /uj
Repeat j
Repeat i
Return
End
4.1.3
RungeKutta Methods
We observed earlier that Heuns method is actually a second-order RungeKutta (RK) method.
Here, we will re-derive Heuns method, via the technique used to derive all explicit RK methods.
The basic idea is to obtain a procedure that agrees with a given Taylor method through a prescribed
order by introducing intermediate function evaluations to replace derivative evaluations. To make
this clear we need to first briefly consider the Taylor series approach.
As always, we begin with the IVP
u = f (u, t) ,
u(t0 ) = u0 .
(4.30)
hn
hn+1
h2
+ + u(n)
+ u(n+1) ()
,
2
n!
(n + 1)!
(4.31)
91
where [t, t + h]. Now, at first glance, this does not appear to be very useful since u is the
unknown function in (4.30). So in particular, it does not appear that derivatives of u should be
known. On the other hand, we have
u (t) = f (u(t), t)
directly from (4.30). Clearly, this can be evaluated at t = t0 ; so the first two terms on the righthand side of (4.31) are known. We should immediately recognize that if we truncate the series after
the second term we will have obtained Eulers method, a (locally) second-order RK scheme.
u ( t)
u( t + h )
u ( t + _2h )
u (t )
t + _2h
t+h
un+1 = un + un h + un
h2
.
2
(4.33)
92
un =
Hence, we can express (4.33) as
un+1
h2
= un + hf (un , tn ) +
2
un + h(un , tn , h) .
f
f
f+
u
t
(un , tn )
(4.34)
(We note in passing that every single-step method can be expressed in the form (4.34).) Our goal
now is to construct an approximation to that differs from in terms of O(h2 ) and higher, and
which does not involve u .
We begin by expanding f in a formal Taylor series as a function of two independent variables.
Thus, we have
f
f
u +
t
f (u + u, t + t) = f (u, t) +
u
t
2
1 f
2f
2f
2
2
+
(u) + 2
ut + 2 (t) .
2 u2
u t
t
(4.35)
f
f
(un , tn ) + ph (un , tn ) + O(h2 ) . (4.36)
u
t
(4.37)
We see that this differs from (4.37) by terms of order h2 , and (4.37) depends only on the function f
and some as yet undetermined constants. We find these constants by comparing (4.38) and (4.34).
Hence,
1
1
a1 + a2 = 1 , a2 p = , and a2 q = .
2
2
93
From the last two of these we see that p = q = 2a12 , provided a2 6= 0, and the first relation yields
a1 = 1 a2 . Thus, given a2 6= 0, we can calculate p, q and a1 , and these constants determine a
function which differs from by terms of order h2 , and higher. We see from this that Runge
Kutta methods are not uniquely determined simply by specifying the order. This can be used to
great advantage by selecting coefficients that, e.g., minimize truncation error.
There are two choices of a2 in wide use. The first is a2 = 21 . Then a1 = 12 , and p = q = 1; so
in (4.37) becomes
(un , tn , h) =
1
1
f (un , tn ) + f (un + hf (un + hf (un , tn ), tn + h) .
2
2
Thus (4.34) is
un+1 = un +
h
[f (un , tn ) + f (un + hf (un , tn ), tn + h)] ,
2
(4.39)
h
h
un + f (un , tn ), tn +
2
2
(4.40)
which is known variously as modified Euler, modified Euler-Cauchy, or simply the midpoint method.
It is of interest to compare the truncation error for these two methods. From (4.35) we see that
the O(h2 ) term for a general RK method has the form
2
2f
2f
h2 2
2 f
(, ) + 2pqf (, )
(, ) + p
(, ) ,
q f (, )
a2
2
u2
ut
t2
with , [u, u + u] [t, t + t]. For Heuns method p = q = 1, and a2 = 21 ; so the truncation
error term is 12 of the form given above. For the midpoint formulas p = q = 21 , and a2 = 1, which
leads to a factor of 41 . Thus, the leading truncation error of the midpoint rule is only half of that
of Heuns method.
We close this section on RungeKutta methods by noting that the higher-order RK schemes
are derived by the same procedure as presented here. But this becomes extremely tedious for
methods of order higher than two. In general, all RungeKutta methods can be written in the form
un+1 = un + h(un , tn , f, h),
where
=
M
X
wi ki ,
i=1
with
ki = f un +
i1
X
j=1
aij kj , tn + ai h .
94
AL
wT
where c is a column vector with first component equal to zero, wT is a row vector, and AL is a
lower triangular matrix. The general representation for a fourth-order RK method is then
0
c2
c3
c4
a21
a31
a41
w1
a32
a42
w2
a43
w3
w4
1
2
1
2
0
0
1
6
1
3
1
3
1
6
Tables of this form can be found for a great variety of RK methods in Lapidus and Seinfield [19].
They provide a useful way to store RK coefficients in a general RK algorithm implementation.
4.1.4
One of the major disadvantages of the RK methods is that as the order of the method is increased
the required number of function evaluations also increases. Hence, if f (u, t) is a complicated vector
function, required arithmetic increases rapidly with order of the method. At least a partial remedy
to this difficulty can be gained by using the multi-step methods to be discussed in this section.
There are many such methods, and many different ways by which they may be derived. Here, we
will restrict attention to the class of explicit procedures known as the AdamsBashforth methods
and the related class of implicit schemes known as the AdamsMoulton methods. These are often
used in combination to construct predictor-corrector techniques. Our approach to deriving these
methods will be similar to that used for the second-order RungeKutta method of the previous
section, namely via the method of undetermined coefficients. However, we will begin with an
outline of the general technique which would be used in a more rigorous construction of these
methods.
As the name implies, a multi-step method is advanced in time using information from more
than just a single previous step. In general, we can consider k-step methods to solve the by now
familiar IVP
u = f (u, t), u(t0 ) = u0 .
If we assume a constant step size h, then formal integration of this equation between time tn and
tn+1 yields
Z tn+1
Z tn+1
f (u, t) dt ,
u dt =
tn
tn
95
tn+1
f (u, t) dt .
(4.41)
tn
Up to this point everything is exactly the same as in our earlier derivation of the trapezoidal
rule. (In fact, the trapezoidal rule is an implicit multi-step methodbut a somewhat trivial one.)
We can evaluate the integral on the right-hand side by replacing f with an interpolation polynomial
spanning the time steps from n k + 1 to n + 1 to obtain a (generally implicit) k-step method.
In particular, we obtain from this a quadrature scheme for approximating the integral in (4.41); so
from our earlier discussion of quadrature methods we expect (4.41) to take the form
un+1 = un + h
k
X
i fn+1i .
(4.42)
i=0
k
X
i un+1i + h
k
X
i fn+1i ,
(4.43)
i=0
i=1
but for our purposes here we will always have 1 = 1, and i = 0, i = 2, . . . , k; so (4.42) will be
the only form of multi-step methods to be constructed herein. Furthermore, we note that whenever
0 = 0, (4.42) is explicit, but otherwise it is implicit and must be solved for un+1 by iteration. We
will now derive i s for both an explicit and an implicit method.
AdamsBashforth Method
For the explicit method we can write (4.42) as
un+1 = un + h
k
X
i fn+1i .
i=1
(4.44)
Since there are actually three coefficients to be determined (including 1 ) we would expect to
produce a method that is exact for polynomials of degree two, and is thus accurate to third order.
Thus, we expand un in a Taylor series to third order about un+1 . (Note that we could also expand
un+1 about un and obtain the same results). This yields
un = un+1 un+1 h + un+1
h2
+ O(h3 ) .
2
Furthermore, we have
un = un+1 un+1 h + u
n+1
h2
+ O(h3 ) = fn ,
2
and
2
3
un1 = un+1 2un+1 h + 2u
n+1 h + O(h ) = fn1 .
96
h2
+
2
h 1 (un+1 un+1 h) + 2 (un+1 2un+1 h) + O(h3 ) .
Finally, we equate powers of h on the right- and left-hand sides to obtain the two equations for
determination of 1 and 2 :
1 + 2 = 1
1
1 + 22 = .
2
The solution to this system is 1 = 32 , 2 = 21 . These are the coefficients of the second-order
(global order is order of the method) AdamsBashforth formula. The higher-order methods can be
derived in exactly the same way.
AdamsMoulton Method
We next carry out the same procedure for an implicit k-step AdamsMoulton method, again for
k = 2. In this case (4.42) is
un+1 = un + h (0 fn+1 + 1 fn + 2 fn1 ) .
(4.45)
We see now that for an implicit k-step method there are four (including 1 ) coefficients to be
determined. As a consequence we obtain an increase of one power of h in the formal local accuracy.
Thus, the Taylor expansions employed must be carried one term further. We have, for example,
un = un+1 un+1 h + un+1
h3
h2
u
+ O(h4 ) .
n+1
2
6
The remaining two expansions given earlier are sufficient. Substitution into (4.45) then yields
h2
h3
un+1 = un+1 un+1 h + un+1 u
n+1
2
6
2
h
+ 2 un+1 2un+1 h + 2un+1 h
+ h 0 un+1 + 1 un+1 un+1 h + un+1
2
1
1 22 un+1 h2
= un+1 + (0 + 1 2 1) un+1 h +
2
1
1
3
4
+
+ 22
u
n+1 h + O(h ) .
2 1
6
The system of equations for the s is then
0 + 1 + 2 = 1
1
1 + 22 =
2
1
1 + 42 = ,
3
97
and we find
2
1
5
, 1 = , 2 =
.
12
3
12
Both AdamsBashforth (4.44) and AdamsMoulton (4.45) methods can be used separately. But
it is typical to use an AdamsBashforth method as an explicit predictor to begin iterations of an
AdamsMoulton method. One can use various combinations of order and k for such a predictorcorrector pair, but it is quite often that a k-step method, with k the same for both predictor
and corrector is used. This is because the order of the complete method is always the order of
the corrector, provided it is iterated to convergence. Moreover, in any case, each iteration of the
corrector results in a gain of one order in h until the order of the corrector is achieved. Since at
least one iteration of the corrector will always be used, it is reasonable to use a predictor that is
one order less accurate and hence, also a k-step method with k the same as the corrector.
0 =
Predictor-Corrector Algorithm
We now provide an algorithm for implementing an arbitrary k-step predictor-corrector pair of
AdamsBashforth/AdamsMoulton type. There are three specific implementational details to be
noted. The first occurs in step 2 where we require generation of k 1 starting values using a
single-step method. The need to do this is one of the major shortcomings of multi-step methods,
in general. In practice, if fourth-order Adams multi-step methods are being used, then fourthorder RungeKutta is typically employed to compute the starting values. But other approaches
are possible.
The second thing to note is that there is really only one term of the implicit method that must
be re-evaluated at each iteration. Thus, a significant amount of arithmetic can be saved if at the
start of each new time step the remaining terms are calculated and placed in a temporary array
to be re-used at each iteration. The elements of this array are denoted gi , i = 1, . . . , neqns in the
algorithm.
Finally, we observe that only a simple fixed-point iteration has been employed for solution of
the implicit difference equations, rather than Newtons method, as was used for the trapezoidal
rule. This is usual for the AdamsMoulton method because, in general, for methods of order higher
than two, the stability properties are not especially good. Thus, the class of problems to which
these methods should be applied is somewhat restricted, and for such problems class simple Picard
iteration converges very rapidly, usually within a few iterations.
Algorithm 4.2 (AdamsBashforth/AdamsMoulton Integrator)
1. Read control flags: neqns, maxitr, nsteps, kAB , kAM
Read numerical parameters: h,
Read initial data: (ui,0 , i = 1, . . . , neqns), t0
k = max (kAB , kAM )
2. Generate k 1 starting values for each of the i = 1, . . . , neqns equations using a single-step
method; e.g., a kth -order RungeKutta method.
3. Begin time stepping
tk1 = t0 + (k 1)h
Do n = k 1, nsteps
4. Evaluate the explicit AdamsBashforth Predictor Formula
Do i = 1, neqns
98
ui,n+1 = ui,n+1
Repeat i
Do m = 1, maxitr
umax = 0
Do i = 1, neqns
(m)
fi,n+1 = f i, ui,n+1 , tn+1
(m+1)
ui,n+1 = gi + h0 fi,n+1
(m+1)
(m)
If ui,n+1 ui,n+1 > umax then
Repeat i
(m+1)
(m)
umax = ui,n+1 ui,n+1
End
iteration error =
99
In closing this section we note that stability analysis for multi-step methods is somewhat more
complicated than for single-step methods. In particular, single-step methods can be adequately
characterized by their region of absolute stability in the h-plane. However, for multi-step methods
it is not sufficient to consider only absolute stability because the difference equations for multi-step
schemes possess multiple solutions. It is always the dominant one that corresponds to the solution
of the differential equation; but in long time integrations, unless the difference approximation is
such that the parasitic solutions remain small compared with the dominant one, the error of the
numerical solution will eventually become unacceptable. It is growth of the parasitic solution that
makes use of multi-step methods difficult for the class of so-called stiff equations which we treat in
the next section.
4.1.5
One of the more difficult classes of problems in numerical IVPs is solution of stiff equations and
stiff systems. These problems arise from a variety of physical situations, but probably were first
identified in chemical kinetics. They are characterized, intuitively, as having at least two widely
disparate time scales. But there is a precise mathematical definition that specifically excludes
problems that are actually unstable. Often this distinction is not made; but it is important because
stability of the underlying differential problem is necessary to the success of any numerical scheme
if long-time integrations are attempted.
The main feature resulting from widely different time scales is that one component of the
solution may quite rapidly decay to zero, and be completely unimportant, while other components
may be slowly varying. One would expect to be able to employ large time steps, and still be able
to predict the slow components with good accuracy. But it turns out that stability is determined
by the rapidly changing component(s); so even though they may be too small to be of importance,
they may nevertheless impose the requirement of very small time steps to maintain stability.
It should also be noted that stiffness is not always so simple as the above description might
imply. In particular, it is possible for all components of a solution to contribute to stiffness, even
though the qualitative behavior of none of them indicates rapid decay. An example of this is given
by the following problem, which can be found in Gear [10], and elsewhere. Consider the linear
system
u = 998u + 1998v
u(0) = 1
v = 999u 1999v
v(0) = 0 ,
v(t) = et + e1000t .
Notice that both components have terms that very rapidly decay to zero, but the solution itself,
changes only slowly, as shown in Fig. 4.5. On the other hand |u | and |v | are large and v 0
at early times. This is what causes the stiffness of this system. This example also indicates that,
contrary to an often-held belief, stiffness does not necessarily arise as a consequence of nonlinearity;
this system is linear.
Formally, for an IVP to be considered stiff, it must have the properties given in the following
definition.
Definition 4.2 The ordinary differential equation initial-value problem for u = f (u, t), u, f Rp ,
is said to be (locally) stiff whenever any eigenvalue of Ju (f ) has Re 0.
100
2
u(t)
1
0
v(t)
1
0.0
0.01
0.02
0.03
0.04
0.05
Observe that this implies u = u, < 0, and || 1 is, in a sense, stiff; that is, stiffness
can occur even for single equations with a single time scale. Recall that this was precisely the
problem with which we demonstrated the difficulties encountered with Eulers method, and the
remedy of the implicit backward Euler method. It can be shown that the trapezoidal integration
exhibits stability properties similar to those of backward Euler, and as we have already seen, it is
considerably more accurate. In particular, both of these methods have the property of A-stability,
defined below, which is theoretically very desirable (but often not very practical) for methods
employed for solving stiff equations.
Definition 4.3 A method is said to be A-stable if all numerical approximations tend to zero as the
number of time steps n when it is applied to the differential equation u = u with a fixed
step size h > 0, and a (complex) constant with Re < 0.
Because the only restrictions on h and are those stated, the implication of the definition is
that a method has an absolute stability boundary at in the left-half complex h-plane if it is
A-stable. The stability diagram for backward Euler is shown in Fig. 4.6, and indicates that this is
true. This can be easily derived from the form of the amplification factor of backward Euler given
earlier in Eq. (4.22).
Similarly, it can be readily checked that the amplification factor for the trapezoidal rule in this
case is
1 + 1 h
2
.
1 12 h
For stability we require that this be less than or equal to unity; thus,
1 + 1 h 1 1 h
2
2
must hold. Clearly, if h > 0 and R with 0, this always holds. If C we have = + i,
101
1
1
1
1
1 + h + i h 1 h i h
2
2
2
2
Stable
Unstable
(1,0)
Stable
Re h
Stable
4.2
In this section we will consider the solution of boundary-value problems (BVPs) for ODEs. These
arise in the analysis of many different problems in engineering and mathematical physics. In the
102
great majority of cases the problems involve a second-order differential equation with boundary
conditions prescribed at each of two ends of a finite interval [a, b]. This is often called a two-point
boundary-value problem. This is the only case we shall treat here; the reader is referred to more
advanced treatments, such as Keller [17] for discussion of other types of BVPs.
There are two widely-used general classes of procedures for solving ODE BVPs. These are i)
shooting, and ii) finite-difference methods. We shall discuss the first of these briefly and present the
second in considerable detail. We then conclude with a somewhat different approach, the Galerkin
procedure, that is more often applied to PDEs but is most easily understood in the ODE context.
We note that there are many other somewhat less-often used approaches, treatments of which are
omitted here for the sake of brevity.
4.2.1
Mathematical Background
We begin by considering some elementary mathematical background material for ODE BVPs. The
most general form of the problem to be considered here is
Lu = f (x, u, u ),
x [a, b] ,
(4.46)
(4.47a)
Bb u(b) = B
(4.47b)
where Ba and Bb are (at most) first-order, but possibly nonlinear, operators. As one might infer
from our earlier discussion of coordinate transformations employed with Gauss quadrature, there
is no loss in generality in taking a = 0 and b = 1, and we will sometimes employ this slight
simplification. We view f as a generally nonlinear function of u and u , but for the present, we will
take f = f (x) only. L is assumed to be a linear second-order operator, so for this case (4.46) is
linear and becomes
Lu a2 (x)u + a1 (x)u + a0 (x)u = f (x) .
(4.48)
In our discussions here, the boundary operators in (4.47) will be such as to lead to one of the
following three types of conditions, here applied at x = a:
Bu(a) = u(a) = A ,
Bu(a) = u (a) = A ,
< 0.
(Dirichlet)
(4.49a)
(Neumann)
(4.49b)
(Robin)
(4.49c)
The same types also apply at x = b where > 0 must hold for the Robin condition. It should
be observed that more general boundary conditions sometimes occur, for example, periodicity
conditions. The reader is referred to Keller [17] for methods of treatment. Also, we note that the
conditions presented in Eqs. (4.49) are often termed first kind, second kind and third kind,
respectively.
A complete boundary-value problem consists of an equation of the form (4.46) (or (4.48)) and
a boundary condition of the form of one of (4.49) applied at each end of the interval. We note that
a different type of boundary condition may be applied at each end, and that in general the value
of A is different at the two ends. For such a problem to have a solution it is generally necessary
either that f (x) 6 0 hold, or that A 6= 0 at one or both ends of the interval. When f (x) 0, and
A = 0 at both ends of [a, b] the BVP is said to be homogeneous and will in general have only the
trivial solution, u(x) 0.
103
B 0 u = B1 u = 0 ,
x [a, b] .
For specific values of (the eigenvalues) this problem will have nontrivial solutions (the eigenfunctions) analogous to the eigenvalues and eigenvectors of linear algebraic systems studied in Chap.
1. We shall not specifically treat eigenvalue problems for differential operators here except to note
that use of finite-difference approximation, as discussed below, reduces this problem to an algebraic
eigenvalue problem to which the methods of Chap. 1 may be applied.
4.2.2
Shooting Methods
The earliest numerical procedures applied to solution of two-point boundary-value problems were
methods based on initial-value problem techniques, such as discussed earlier. The basic notion is
to view (4.48) as an initial-value problem starting at x = 0. Since (4.48) is second order, two
initial conditions are required. But only one boundary condition can be given at x = 0. Suppose
it is the Dirichlet condition u(0) = A. We must now guess a value for u (0) in order to start the
initial-value problem. Take u (0) = (1) , and integrate (4.48) from x = 0 to x = 1 using any
stable initial-value method. The result might appear as in Fig. 4.7. In particular, the boundary
condition, say u(1) = B, will not have been satisfied. So we select a second value of u (0), say (2) ,
and shoot again. The original curve suggests that we should take (2) > (1) ; the result might
be as shown.
This approach can be formalized as a Newton iteration by observing that the value of u(1) that
is actually obtained is an implicit function of . In particular, we can employ Newtons method to
solve the equation
F () = u(1, ) B = 0 ,
for . For linear equations this turns to out to work fairly well because an auxiliary equation can
be derived in such a way that exactly two iterations are always required (see [17]). However, the
situation is not nearly so favorable for nonlinear problems.
(2)
u= B
(2)
u= B
u= B
(1)
u= A
x=0
(1)
x=1
104
In general, shooting is felt to have the advantage of high-order accuracy available via the welldeveloped initial-value methods. But it has the disadvantage of being iterative, thus requiring
implementation of some algebraic equation solving method in addition to the IVP method. In
addition, the method may completely fail for problems whose solutions exhibit steep gradients
anywhere in the domain.
4.2.3
Finite-Difference Methods
For the above reasons, finite-difference methods have become more widely used in recent years.
It was at one time felt that the standard centered-difference approximations to be employed here
were often not sufficiently accurate. However, when these are used in conjunction with Richardson
extrapolation, the resulting solutions are accurate to fourth order, which is equivalent to the accuracy attained with usual implementations of the shooting method (typically using a RungeKutta
initial-value solver). Furthermore, the finite-difference methods still work satisfactorily in regions
of large gradients, at least if local mesh refinement is employed.
The basic idea in applying the finite-difference method is extremely simple. We merely replace
all of the differential operators in (4.48) and (4.49) with difference operators as treated in Chap. 3.
This results in a system of algebraic equations which, as we will see below, is banded and sparse
(in fact, tridiagonal). Thus, it can be very efficiently solved by the methods of Chap. 1.
Discretization of the Differential Equation
We begin by introducing a partition of the interval [0, 1] consisting of N points: {xi }N
i=1 , such that
0 = x1 < x2 < < xN 1 < xN = 1. Furthermore, we shall always employ a uniform spacing
for the computational mesh. When local mesh refinement is needed, it is best accomplished via a
coordinate-stretching transformation prior to discretization of the domain and equations. However,
we shall not treat this matter here. The interested reader is referred to Thompson et al. [35]. Figure
4.8 depicts the computational mesh. The solution to (4.48) will be approximated at each of the
N grid points, and the ordered collection of these grid point values, {ui }N
i=1 , constitutes the grid
function. Finally, we observe that for the general interval [a, b], we would have h = (b a)/(N 1).
We will later see that grid points i = 0 and i = N + 1 are needed for some types of boundary
0
1
x=0
3 i1
2
h
i
x = (i1)h
i+1
N1
N+1
x=1
(4.50)
One expects, correctly, that this approximation exhibits a second-order truncation error; moreover,
the expansion for the truncation error contains only even powers of h. To demonstrate this, we
first carry out the indicated differencing. We have
ui+1 ui1
ui1 2ui + ui+1
+ a1,i
+ a0,i ui = fi .
a2,i
h2
2h
105
h
a1,i (ui+1 ui1 ) + a0,i h2 ui = h2 fi ,
2
(4.51)
In general, this difference equation holds for all interior grid points on (0, 1); i.e., for i = 2, . . . , N 1.
Some modifications, depending on boundary condition type, will be needed at the two boundary
points. We treat this later.
For notational convenience, we now define
h
,
2
a0,i h2 2a2,i ,
h
a2,i + a1,i .
2
(4.52)
This form strongly suggests the tridiagonal matrix structure alluded to earlier. But it will be
necessary to show that introduction of boundary conditions does not alter this structure.
Truncation Error Analysis
We will next derive the truncation error for the interior grid point approximations (4.52). As the
reader should expect by now, the approach makes use of Taylor expansions, and as a consequence,
a certain degree of smoothness will be required of u(x) in order for the truncation error estimates
to be valid. For ui1 and ui+1 expanded about the arbitrary grid point xi we have
ui1 = ui ui h + ui
h2
h3
h4
u
+ u
,
i
i
2
6
24
and
4
h3
h2
h
+ u
+
u
+ .
i
i
2
6
24
Substitution of these into (4.52), and rearrangement by grouping terms containing like powers of
h, leads to
ui+1 = ui + ui h + ui
h2
u
2 i
h4
h3
+
(C
+
C
)
u + = h2 fi .
+ (C3,i C1,i ) u
1,i
3,i
6 i
24 i
106
h4
h4
ui + a2,i u
= h2 fi .
6
12 i
107
(4.55)
The main observation to be made here is that the tridiagonal matrix structure is now maintained,
and at the same time second-order accuracy will be achieved. It is left to the reader to derive a
similar result at i = N , where an image point corresponding to i = N + 1 occurs. Finally we note
that the same image points occur for Robin boundary conditions, and that they may be treated in
precisely the same manner as shown here for the Neumann condition. We leave demonstration of
this as an exercise for the reader.
Complete BVP Discretization
We summarize the treatment given here by presenting the matrix structure and pseudo-language
algorithm for solving the following problem.
a2 (x)u + a1 (x)u + a0 (x)u = f (x),
x [0, 1) ,
u (0) = A ,
(4.56)
u(1) = B .
The first equation of the system will be (4.55). Then the next N 2 equations will all be of the
form (4.52), and the last equation will be of the form (4.53) but written for the grid point xN .
The matrix form of this is shown in Fig. 4.9. The tridiagonal structure is quite evident from this
u1
C2,1 C1,1 +C3,1
f1 + h2 AC1,1
u2
C1,2
f2
C2,2
C3,2
f3
C1,3
C2,3 C3,3
u3
..
..
..
..
..
.
.
.
.
. = h2
fi
C1,i C2,i C3,i
ui
..
..
..
..
..
.
.
.
.
.
fN 1
C1,N 1 C2,N 1 C3,N 1 uN 1
B/h2
0
1
uN
Figure 4.9: Matrix structure of discrete equations approximating (4.56)
figure. Clearly, the system can be efficiently solved for the grid function values {ui }N
i=1 using the
tridiagonal LU-decomposition algorithm of Chap. 1. A pseudo-language algorithm for implementing
this approach follows.
Algorithm 4.3 (Finite-Difference Solution of Linear BVPs)
1. Enter problem data: Number of grid points, N ; endpoints of solution domain, a, b; boundary
condition type flags, ibca, ibcb; boundary values A, B.
2. Calculate discretization step size: h = (b a)/(N 1).
3. Load matrix coefficients for all grid points
Do i = 1, N
108
4.3
In this section we will consider several special topics, any or all of which can be of great utility in
applications. We begin with a brief discussion of coordinate singularities, as arise from use of polar
coordinate systems. This is followed by a basic treatment of nonlinear equations. Four approaches
are presented: i) Picard iteration, ii) Newtons method, iii) quasilinearization, and iv) Galerkin
procedures.
4.3.1
109
Coordinate Singularities
Correct treatment of coordinate singularities is important for problems posed in polar geometries;
consider the following model problem:
1
u + u u = f (r) ,
r
u (0) = 0 ,
r [0, R) ,
u(R) = 1 .
We begin by noting that this is in the same general form as the boundary-value problems already
treated, and thus, away from singularities approximations will be the same as those of the preceding
section; but clearly, when r = 0, the coefficient of the first-derivative term is singular. We might
first consider multiplying through by r, but this simply yields the boundary condition at r = 0,
and we have no guarantee that the differential equation, itself, will be satisfied. In this particular
situation, it is typical that the condition u (0) = 0 results from geometric symmetries, and is not
even an actual boundary condition since r = 0 is, in fact, not really a boundary point in such a
case. Thus, it is crucial to guarantee that the differential equation also be satisfied.
This requirement can be met by considering the limit of the differential equation as r 0. We
have
1
u 0 = u2 .
110
Hence, the tridiagonal form will still be preserved while satisfying both the differential equation
and symmetry condition to second order.
We should remark here that although the above development is straightforward and mathematically consistent, other approaches have been widely used; these often yield reasonable results,
but can sometimes completely fail unexpectedlyhence, we do not recommend them. One such
approach is to implement only the Neumann boundary condition, and ignore the differential equation, at the point of singularity, as mentioned earlier. Now suppose in our preceding example we
imposed the condition u(R) = 0, and in addition defined f (r) as
(
1 r = 0,
f (r) =
0 otherwise .
This would correspond to a steady-state heat conduction problem in a cylinder of radius R with
heating due to a line source on the axis. If we were to implement only the Neumann condition at
r = 0, thus avoiding the singularity, we would no longer have any heat input, and the solution would
be the trivial one, u 0. The only approach that can consistently circumvent such difficulties is
the one utilizing LHospitals rule for removal of the coordinate singularity, as we have described
above.
4.3.2
Our next topic is solution of nonlinear two-point boundary-value problems. Recall that our original
equation (4.46) was formally nonlinear. We repeat it here:
Lu = f (x, u, u ) ,
(4.57)
where, again, L is a linear operator of order two. We will begin by considering three methods
for solving (4.57) in the context of finite-difference approximation: i) Picard iteration, ii) Newton
iteration, and iii) quasilinearization. We then conclude our discussions of ODE BVPs in the next
section by presenting a Galerkin approximation for a specific case of (4.57).
Picard Iteration
If we replace L with Lh , where from earlier notation
Lh a2,i D02 (h) + a1,i D0 (h) + a0,i ,
then, as we have already seen, the left-hand side of the resulting system of discrete equations is
just a matrix times a vector, as shown in detail following Eq. (4.56). Thus, we can write the formal
solution representation of (4.57) as
u = L1
h f (x, u, u ) ,
(4.58)
111
We know that a sufficient condition for convergence of this iteration procedure is that the
Lipschitz constant be less than unity. For many mild nonlinearities in f , this can be proven to
hold. In such cases, the solution of (4.57) is straightforward. We simply discretize the equation as
ui+1 ui1
,
Lh ui = f xi , ui ,
2h
i = 1, 2, . . . , N ,
(4.59)
ui+1 ui1
xi , ui ,
2h
=0.
(4.60)
As the notation implies, each of the Fi s depends on only three grid function values because of our
use of second-order centered differences. It then follows that J(F ), the Jacobian matrix of F , is
tridiagonal; so the Gaussian elimination procedure usually employed for Newtons method should
be replaced with the tridiagonal LU-decomposition algorithm of Chap. 1. We will demonstrate
these ideas with the following example.
Consider the problem
u = S(x) (u )2 sin u f (x, u, u ) ,
x (0, 1) ,
u(0) = u(1) = 0 .
(4.61)
Here, S(x) is a known forcing function. The differential equation then has L = d2 /dx2 , and the
second-order accurate centered-difference approximation is
ui1 2ui + ui+1
= Si
h2
ui+1 ui1
2h
2
sin ui .
As we have indicated earlier in the linear case, it is preferable to multiply through by h2 , so this
becomes
1
ui1 2ui + ui+1 = h2 Si (ui+1 ui1 )2 h2 sin ui .
4
Then for each i = 2, . . . , N 1, we have
1
Fi = ui1 2ui + ui+1 + (ui+1 ui1 )2 + h2 sin ui h2 Si = 0 ,
4
(4.62)
(4.63)
112
The solution procedure for the nonlinear BVP (4.61) now consists of solving the system of
nonlinear algebraic equations, (4.62), (4.63) via Newtons method. We compute the elements of
the Jacobian matrix of F in the usual way. For general i, we have
Fi
u1
Fi
ui1
Fi
ui
Fi
ui+1
Fi
ui+2
=0,
Fi
=0,
u2
... ,
Fi
=0,
ui2
1
= 1 (ui+1 ui1 ) ,
2
= 2 + h2 cos ui ,
1
= 1 + (ui+1 ui1 ) ,
2
Fi
= 0 , ... ,
=0,
uN
i = 2, . . . , N 1 .
For i = 1, F1 /u1 = 1, and all other partial derivatives are zero; similarly, FN /uN = 1, with all
other partial derivatives being zero. It is clear that the Jacobian matrix is tridiagonal as we noted
earlier. Thus, Newtons method can be implemented in a very efficient manner.
Quasilinearization
Quasilinearization is generally equivalent to Newtons method, except in the specifics of its implementation. In Newtons method, as just presented, we discretize the nonlinear equation and
then locally linearize the nonlinear (algebraic) difference equations at each iteration. By way of
contrast, in quasilinearization (often called the NewtonKantorovich procedure) we first linearize
the nonlinear operator(s), and then discretize the resulting linear equation(s). Iterations are then
performed in a manner similar to that of Picard iteration, but of course, using a different iteration
function that results in quadratic convergence just as for Newtons method.
Recall that our prototypical nonlinear equation is (4.57):
Lu = f (x, u, u ) ,
where f is generally nonlinear in u and u . We view u and u as being distinct independent functions
and expand f in a Taylor series, called a FrechetTaylor series, in terms of these:
f (0)
f (0)
(0)
(0)
(0) (0)
+ ,
+
u
u
uu
+
f (x, u, u ) = f x, u , u
u
u
where u(0) and u (0) are initial estimates of u and u , respectively, which are typically updated after
each iteration. Equation (4.57) now becomes
"
(0) #
f (0)
f (0) (0)
f
f (0) d
(0) (0)
(0)
u
.
(4.64)
u = f x, u , u
u
L
u
dx
u
u
u
It is clear that this equation is linear in u, and that it can be discretized in the usual manner.
In particular, the left-hand side can still be cast in the form of (4.52). To see this we write the
left-hand side in detail:
"
"
(0)
(0) #
(0)
(0) #
f
f
f
f
u
u + a0
u = a2 u + a1
u.
a2 u + a1 u + a0 u
u
u
u
u
113
u
,
u
f(x) f x, u , u
u
u
we obtain a linear equation of precisely the same form as treated earlier; namely
a
2 (x)u + a
1 (x)u + a
0 (x)u = f(x) .
(4.65)
As a consequence of the preceding definitions, the difference equation coefficients at the mth iteration
take the form
"
(m) #
f
h
(m)
C1,i = a2,i a1,i
u i
2
"
(m) #
f
(m)
h2 2a2,i
C2,i = a0,i
u i
"
#
f (m) h
(m)
C3,i = a2,i + a1,i
u i
2
f (m)
f (m)
(m)
(m)
(m)
(m)
(m)
fi = f xi , ui , D0 ui
D0 ui .
ui
u i
u i
In our usual notation for iteration procedures, the discrete form of (4.65) now can be expressed as
(m) (m+1)
C1,i ui1
(m) (m+1)
+ C2,i ui
(m) (m+1)
+ C3,i ui+1
(m)
= h2 fi
at the general ith interior grid point. Clearly, boundary condition implementation is done in
the same way as demonstrated earlier for linear equations. Furthermore, we note that nonlinear
boundary conditions can be linearized via quasilinearization in a manner analogous to the treatment
presented here for the differential equation. The algorithm required to implement quasilinearization
for solution of (4.57) is the following.
Algorithm 4.4 (Quasilinearization Solution of Nonlinear ODE Two-Point BVPs)
1. Input number of grid points N , maxitr, , boundary points a, b, boundary condition flags and
o
n
(0) N
; set iteration counter, m = 0; calculate h = (b a)/(N 1)
values, initial guess ui
i=1
(m)
, ui (m)
i = 1, . . . , N .
(m)
(m)
(m)
(m)
4. Calculate C1,i , C2,i , C3,i , fi , i = 1, 2, . . . , N , and store in banded matrix
114
o
(m+1) N
ui
i=1
We note in closing this section that this function space version of Newtons method exhibits the
same quadratic convergence as does the usual Newtons method on RN .
4.3.3
The Galerkin procedure possesses features that are fundamentally different from finite differencing.
In finite-difference methods, all approximations are local, extending over only a few grid points,
and their derivations are based on Taylor series expansion, thus implying considerable smoothness
requirements on the quantities being approximated. In the Galerkin procedure we begin by assuming a global functional representation of the solution. For reasons that will become more clear as
we proceed, this global representation should be in the form of a (generalized) Fourier series, say
u(x) =
X
k=1
ak k (x) ,
x [a, b] .
(4.66)
The functions {k }
k=1 are called the basis functions of the Galerkin approximation. They are
usually chosen to automatically satisfy any boundary conditions that are to be satisfied by u(x),
although this is not strictly necessary and certainly not always possible.
A more important, and in fact, necessary, requirement is that the set of functions {k (x)}
k=1
be complete in L2 [a, b]. In general, completeness is not an easy property to check. However, it is
known that the set of eigenfunctions of any regular (or even singular) SturmLiouville problem is
complete (see, e.g., Gustafson [11]), and this supplies us with a quite rich assortment of possible
basis sets. The importance of completeness cannot be overstated, for completeness of a given basis
set in L2 [a, b] implies that any function in L2 [a, b] can be represented by a series of the form (4.66).
Hence if the solution to the BVP under consideration is in L2 we are guaranteed that Fourier
coefficients {ak }
k=1 exist, and that the resulting series (4.66) is convergent. Without completeness
of the set {k }, we have no such guarantee.
We observe that there is one other property of {k } that is desirable. It is orthonormality.
Clearly, if the set {k } is complete it can be made orthonormal via the GramSchmidt procedure (see,
e.g., Stakgold [29]) followed by normalization. Orthonormality implies the following relationship
amongst the k s:
hi , j i = ij , i, j = 1, . . . , ,
(4.67)
R
where hi , j i = i j , and ij is the Kronecker . Although this is not an absolute requirement,
it is a useful property in terms of implementing the Galerkin procedure, and also for testing convergence of the series (4.66). It is recommended that orthonormal basis sets essentially always be
used.
Finally, we note that the summation indexing, k = 1, . . . , , is not the only possible one.
Choice of indexing is dictated by details of the problem being solved and the basis set employed to
do this, and could also be either k = 0, . . . , or k = inf ty, . . . , 0, . . . , in place of the one used
herein.
115
Example Problem
It is probably best to illustrate implementation of the Galerkin procedure by means of a simple,
but fairly general example. We choose the following relatively simple, but nonlinear, problem:
u u2 = f (x) ,
x (0, 1) ,
(4.68)
(4.69)
Observe that because this is a nonlinear problem some minor complications will arise in implementing the Galerkin procedure, thus providing an opportunity to discuss some subtleties.
Our first task is to choose a set of basis functions. A general way to do this is to consider
the eigenvalue problem corresponding to the linear part of the differential operator. Thus, we seek
{k }
k=1 satisfying
= , x (0, 1)
(4.70)
(0) = (1) = 0 .
(4.71)
Since this is an eigenvalue problem associated with a very basic SturmLiouville problem, it is
known that the eigenfunctions are complete in L2 [0, 1]. Furthermore, the principal part of the
operator (that part consisting of highest-order derivative terms) is the same in both (4.68) and
(4.70). Thus, we expect a series of the form (4.66) to provide a solution to the problem (4.68, 4.69)
when k s are nontrivial solutions of (4.70, 4.71).
It is also worth mentioning that in the case where the boundary conditions (4.69) are nonhomogeneous, we might still use as a basis the same functions k obtained from (4.70, 4.71). These
would no longer satisfy the boundary conditions, say
u(0) = A ,
u(1) = B ;
(4.72)
but we can make a slight modification to (4.66) to remedy this. Namely, we replace (4.66) with
u(x) = A(1 x) + Bx +
{k }
k=1 .
ak k (x) .
(4.73)
k=1
provided we assume commutativity of differentiation and series summation in the first term.
We note that the nonlinear term requires some special attention. Namely, it should be rewritten
as
!2
!
!
X
X
X
ai i =
ai i
ai i
i
X
i
ai i
X
j
aj j ,
116
to prevent inadvertent loss of terms from the final product series. We next rearrange these so that
!2
X
XX
ai i
=
ai aj i j .
i
In general, the validity of such a rearrangement, and hence of the above indicated equality, requires
absolute convergence of the individual series involved in the rearrangement. However, for the Fourier
series considered here it can be shown, via the Parseval identity (see, e.g., [29]), that less is required.
Indeed, all that is necessary is convergence in 2 . But this is required for (4.66) to be a solution, in
the first place.
We next rewrite (4.74) as
X
X
ai i (x)
ai aj i j = f (x) .
(4.75)
i
i,j
i,j
u u2 , k = hf, k i ,
is called the weak form of Eq. (4.68). The reason is that the ak s obtained from (4.76) do no
necessarily lead to a function u(x), via (4.66), that satisfies Eq. (4.68), itself; that is,
u u2 = f (x).
The reason for this is quite clear. Equation (4.76) is an integral relation while (4.68) is a differential
equation. It is possible that ak s determined from (4.74) will not lead to a function u(x) that
is smooth enough to permit the differentiation required in (4.68). Recall that we only require
{ak } 2 , which implies only that u(x) L2 .
Now we observe that the system of equations, (4.76), is not suitable for numerical computation
because it is of infinite (countable in this case) dimension. To implement the Galerkin procedure
we must be content to consider approximations
uN (x) =
N
X
ak k (x)
k=1
such that
u(x) = uN (x) +
k=N +1
= uN (x) + R(x) .
ak k (x)
(4.77)
117
Clearly, if the series (4.66) is convergent (which is required for existence of a solution) N (depending on ) such that
ku(x) uN (x)k = kR(x)k < > 0 ,
for some norm, k k (in particular, for the L2 -norm). Thus, for some N we can carry out the
preceding development using (4.77) in place of (4.66) and expect to obtain a reasonably good
approximation to the solution of (4.68, 4.69).
In place of the infinite system (4.76), we now have the N N system
+
+ * N
*N
X
X
ai aj i j , k = hf, k i
(4.78)
ai i , k
i,j=1
i=1
for k = 1, 2, . . . , N . We now carry out the indicated integrations (which in some cases may require
numerical quadrature) to obtain
N
X
i=1
or
N
X
i=1
hi , k i,
ai hi , k i
ai Aik
N
X
N
X
i,j=1
ai aj hi j , k i = hf, k i ,
ai aj Bijk = Fk ,
k = 1, 2, . . . , N ,
(4.79)
i,j=1
with Aik
Bijk hi j , k i. Equation (4.79) is a N N nonlinear system for the ak s; it
is customary to solve it via Newtons method treated in Chap. 2.
Convergence of Galerkin Procedures
It is important to note that the matrices corresponding to the values of the inner products, Aik
and Bijk , are nonsparse. So the values of permissible N are generally much smaller than is the case
for finite-difference methods. On the other hand, it can be shown that under certain smoothness
assumptions on u, and with appropriate types of boundary conditions, the Galerkin procedure (and
spectral methods, in general) are infinite order. This means that if we define the error to be
eN ku(x) uN (x)k
for some norm k k, then eN 0 faster than any finite power of
1
eN o
p<.
Np
1
N.
That is,
There are two separate convergence tests that must be performed on Galerkin solutions to
nonlinear problems. The first is convergence of the Newton iterations used in the solution of
(4.79). This can be done in the usual manner, since as far as these iterations are concerned, only a
finite-dimensional system is involved. The second is convergence of the series representation (4.77).
This is analogous to convergence of grid functions in the case of finite-difference solutions. But now
we monitor convergence with respect to increasing the number of terms in the series expansion.
We will consider three alternative criteria by which to test convergence.
The first is the one that we would always prefer to satisfy because it implies that the solution
is as smooth as required by the differential equation. Namely, we check that
eN max uN u2N f < ,
x[0,1]
118
for all N greater than or equal to some value N0 . This is the strong operator max-norm error or
max-norm residual. It is possible, for reasons alluded to earlier, that convergence in this norm may
not be achieved.
The second convergence test is strong convergence in L2 . This implies that eN < for all N
greater than some specific N0 with eN defined in terms of the L2 norm as
eN
Z
1
0
[uN
u2N
f (x)] dx
21
< .
Z
1
0
[u(x) uN (x)] dx
21
<.
This simply implies convergence of the Fourier representation, but it does not show that the differential equation, itself, has been satisfied. Since, in general, we do not know u(x), we typically
replace the above with
21
Z 1
2
<,
[uN +M (x) uN (x)] dx
0
where M 1. If the k s are orthonormal, it follows from (4.77) and Parsevals identity that we
may simply test
NX
+M
a2k < .
k=N +1
4.3.4
Summary
The initial value problem methods treated in this chapter have included both single-step and
multi-step techniques, but with emphasis on low-order versions of the former. This has been done
for simplicity of presentation and also due to the wide use of such methods in the context of
discretizations of partial differential equations to be studied in the next chapter.
We have also covered a broad range of topics associated with two-point boundary-value problems for ODEs. In particular, shooting, finite-difference and Galerkin methods have all received
attention; and in addition special topics including treatment of coordinate singularities and nonlinearities have been discussed. The emphasis, however, has clearly been on finite-difference methods
because these will be employed in our numerical treatment of partial differential equations to follow.
Chapter 5
5.1
Mathematical Introduction
In the present treatment we will mainly restrict attention to the partial differential equations of
classical mathematical physics, namely the heat, wave and Laplace/Poisson equations, and we
will treat these only with well-known, classical numerical methods. This is done for the sake
of brevity; but it is important that these elementary methods be mastered prior to any attempt
to study more complicated problems and correspondingly more advanced numerical procedures.
Moreover, these relatively simple methods can often be valuable in their own right.
In contrast to the situation for ordinary differential equations, the state of development of
numerical methods for partial differential equations is somewhat less advanced. This is so for a
variety of reasons, not the least of which is the fact that, theoretically, PDEs are more difficult to
treat than are ODEs. In turn, this is true for a number of reasons. Recall that there are only two
basic types of problems associated with ODEs: IVPs and BVPs. Moreover, in each of these cases
the required auxiliary conditions are readily identified, and they can be implemented numerically
in a quite uniform manner across a broad spectrum of individual problems.
In the PDE case there are three separate classes of problems: i) IVPs, usually called Cauchy
problems in the PDE context, ii) BVPs, and iii) initial boundary value problems (IBVPs). Furthermore, within each of these classes, individual problems may be formulated for one or more
types of linear differential operators, namely elliptic, parabolic, or hyperbolic, as well as various
less easily classified nonlinear operators. The generic features of solutions, and to some extent
solution techniques, differ among these problem types. Additionally, because partial differential
equations are often posed on domains of more than one spatial dimension, boundary conditions
must usually be imposed along a curve (or even on a surface), rather than simply at two endpoints
of an interval as in the ODE case. Likewise, initial data must be prescribed for the entire spatial
domain, instead of at a single point as for ODE IVPs. All of these features combine to provide
119
120
a richness in variety of PDE problems that is not present for ODE problems. But this richness
implies difficulties that also are not present in the ODE case. Nevertheless, within the framework
of the problems to be considered in these lectures, the theoryboth analytical and numericalis
complete and well understood.
There are two pieces of mathematics regarding partial differential equations that are extremely
important for numerical analysis. The first of these, classification, is one of the first topics covered in any elementary PDE course, while the second, well posedness, is generally introduced in
more advanced courses. Here, we will briefly discuss each of these before proceeding to numerical
treatment of specific classes of problems.
5.1.1
The three above mentioned types of equations are all examples of linear second-order partial differential equations. The most general form of such equations, when restricted to two independent
variables and constant coefficients, is
auxx + buxy + cuyy + dux + euy + f u = g(x, y) ,
(5.1)
where g is a known forcing function; a, b, c, . . . , are given constants, and subscripts denote partial
differentiation. In the homogeneous case, i.e., g 0, this form is reminiscent of the general
quadratic form from high school analytic geometry:
ax2 + bxy + cy 2 + dx + ey + f = 0 .
(5.2)
Equation (5.2) is said to be elliptic, parabolic or hyperbolic according as b2 4ac is less than, equal
to, or greater than zero. This same classification is employed for (5.1), independent of the nature
of g(x, y). In fact, it is clear that the classification of linear PDEs depends only on the coefficients
of the highest-order derivatives. This grouping of terms,
auxx + buxy + cuyy ,
is called the principal part of the differential operator in (5.1), and this notion can be extended in
a natural way to more complicated operators. Thus, the type of a linear equation is completely
determined by its principal part. It should be mentioned that if the coefficients of (5.1) are permitted to vary with x and y, its type may change from point to point within the solution domain.
This can pose difficulties whose treatment requires methods far beyond the scope of the material
to be presented in these lectures.
We next note that corresponding to each of the three types of equations there is a unique
canonical form to which (5.1) can always be reduced. We shall not present the details of the
transformations needed to achieve these reductions, as they can be found in many standard texts
on elementary PDEs (e.g., Berg and MacGregor [2]). On the other hand, it is important to be
aware of the possibility of simplifying (5.1), since this may also simplify the numerical analysis
required to construct a solution algorithm.
It can be shown when b2 4ac < 0, the elliptic case, that (5.1) collapses to the form
uxx + uyy + Au = g(x, y) ,
(5.3)
121
which is the heat equation, or the diffusion equation; and for the hyperbolic case, b2 4ac > 0,
Eq. (5.1) can always be transformed to
uxx uyy + Bu = g(x, y) ,
(5.5)
where B = 0 or 1. If B = 0, we have the wave equation, and when B = 1 we obtain the linear
KleinGordon equation. The Laplace, heat and wave equations are often termed the equations
of classical mathematical physics, and we will consider methods for solving each of these in the
sections that follow. We remark that the other linear equations just mentioned can be solved by
the same methods within each classification.
5.1.2
Before proceeding to introduce numerical methods for solving each of the three main classes of
problems it is worthwhile to give some consideration to the question of under what circumstances
these equations do, or do not, have solutions. This is a part of the mathematical concept of well
posedness. We begin with a formal definition of this property.
Definition 5.1 A problem consisting of a partial differential equation and boundary and/or initial
conditions is said to be well posed in the sense of Hadamard if it satisfies the following conditions:
i) a solution exists;
ii) the solution is unique;
iii) the solution depends continuously on given data.
The well-posedness property is crucial in solving problems by numerical methods because essentially all numerical algorithms embody the tacit assumption that problems to which they apply
are well posed. Consequently, a method is not likely to work correctly on an ill-posed (i.e., a not
well-posed) problem. The result may be failure to obtain a solution; but a more serious outcome
may be generation of numbers that have no association with reality in any sense. It behooves the
user of numerical methods to understand the mathematics of any considered problem sufficiently to
be aware of the possible difficultiesand symptoms of these difficultiesassociated with problems
that are not well posed.
We close this discussion of well posedness by describing a particular problem that is not well
posed. This is the so-called backward heat equation problem. It arises in geophysical studies
in which it is desired to predict the temperature distribution within the Earth at some earlier
geological time by integrating backward from the (presumed known) temperature distribution of
the present. To demonstrate the difficulties that arise we consider a simple one-dimensional heat
equation.
ut = uxx ,
x (, ) , t [T, 0) ,
with
u(x, 0) = f (x) .
Formally, the exact solution is
u(x, t) =
1
4t
f ()e
(x)2
4t
d ,
(5.6)
the derivation of which (see, e.g., Berg and MacGregor [2]) imposes no specific restrictions on the
sign of t. But we immediately see that if t < 0, u(x, t), if it exists at all, is imaginary (since
122
, the thermal diffusivity, is always greater than zero). In fact, unless f decays to zero faster
than exponentially at , there is no solution because the integral in (5.6) does not exist. It
will turn out that this behavior of heat equation solutions will place restrictions on the form of
difference approximations that can be used to numerically solve it. Later we will present a very
natural difference approximation to the heat equation that completely fails, in part, because of
nonexistence of solutions to the backward heat equation.
5.2
Although in the sequel we will consider only basic finite-difference methods for approximating
solutions of partial differential equations, in the present section we shall provide brief discussions
of several of the most widely-used discretization techniques. We note at the outset that temporal
discretizations are almost always either finite-difference or quadrature based, but many different
methods, and combinations of these, may be employed for spatial approximation. Figure 5.1 depicts
the main features of what are the most well-known classes of methods: i) finite difference, ii), finite
element and iii) spectral.
y
(a)
(b)
typical
mesh stencil
i,j
geometric
"element"
j th horizontal
grid line
i th vertical grid line
nodal point
x
y
(c)
no geometric domain
discretization required
Figure 5.1: Methods for spatial discretization of partial differential equations; (a) finite difference,
(b) finite element and (c) spectral.
As we will describe in considerable detail in later sections, finite-difference methods are con-
123
structed by first gridding the solution domain as indicated in Fig. 5.1(a), and then deriving
systems of algebraic equations for grid-point values which serve as approximations to the true solution at the discrete set of points defined (typically) by intersections of the grid lines. The figure
depicts the 2-D analogue of what was already been done for ODE BVPs in Chap. 4. We remark
that the regular structured grid shown here is not the only possibility. But we will not treat
finite-difference approximations on unstructured grids in the present lectures. (The interested
reader is referred to Thompson et al. [35] for comprehensive treatments of both structured and
unstructured gridding techniques, generally referred to as grid generation.)
Finite-element methods (FEMs) are somewhat similar to finite-difference methods, although
they are typically more elaborate, often somewhat more accurate, and essentially always more difficult to implement. As can be seen from Fig. 5.1(b) the problem domain is subdivided into small
regions, often of triangular (or, in 3-D, tetrahedral) shape. On each of these subregions a polynomial (not unlike those discussed in Chap. 3, but multidimensional) is used to approximate the
solution, and various degrees of smoothness of the approximate solution are achieved by requiring
constructions that guarantee continuity of a prescribed number of derivatives across element boundaries. This approximation plus the subregion on which it applies is the element. It should be
remarked that the mathematics of FEMs is highly developed and is based on variational principles
and weak solutions (see, e.g., Strang and Fix [31]), in contrast to the Taylor-expansion foundations
of finite-difference approaches. In this sense it bears similarities to the Galerkin procedure discussed
in Chap. 4.
Figure 5.1(c) displays the situation for spectral methods. One should observe that in this case
there is no grid or discrete point set. Instead of employing what are essentially local polynomial
approximations as done in finite-difference and finite-element methods, assumed forms of the solution function are constructed as (generalized) Fourier representations that are valid globally over
the whole solution domain. In this case, the unknowns to be calculated consist of a finite number
of Fourier coefficients leading to what amounts to a projection of the true solution onto a finitedimensional space. We previously constructed such a method in Chap. 4 in the solution of ODE
BVPs via Galerkin procedures. In modern terminology, these are examples of so-called grid-free
methods.
There are numerous other discretization methods that have attained fairly wide acceptance in
certain areas. In computational fluid dynamics (CFD) finite-volume techniques are often used.
These can be viewed as lying between finite-difference and finite-element methods in structure,
but in fact they are essentially identical to finite-difference methods despite the rather different
viewpoint (integral formulation of governing equations) employed for their construction. Other less
often-used approaches deserving of at least mention are spectral-element and boundary-element
methods. As their names suggest, these are also related to finite-element methods, and particularly
in the latter case are applicable mainly for a very restricted class of problems. Beyond these are a
number of different pseudo-spectral methods that are similar to spectral methods, but for which
the basis functions employed are not necessarily eigenfunctions of the principal part of the differential operator(s) of the problem being considered. In addition to all these individual techniques,
various calculations have been performed using two or more of these methods in combination.
Particular examples of this include finite-difference/spectral and finite-element/boundary-element
methods. Discussion of details of such schemes is far beyond the intended scope of the present lectures, and from this point onward our attention will be focused on basic finite-difference methods
that have been widely used because of their inherent generality and relative ease of application.
124
5.3
Parabolic Equations
In this section we will consider several widely-used methods for solving parabolic problems. We
begin with forward- and backward-Euler methods applied to the 1-D heat equation, the former
for the Cauchy problem and the latter for an initial boundary value problem. Both of these are
first-order accurate in time. We then treat two second-order methods. The first is unconditionally
unstable and should never be used, but it is a very intuitively natural scheme; the second is the
well-known unconditionally stable CrankNicolson method. We end the section with analysis of the
PeacemanRachford alternating direction implicit method, a very efficient technique for applying
the CrankNicolson scheme to 2-D parabolic equations.
5.3.1
We begin our study of the numerical solution of parabolic equations by introducing the simplest
procedure for solving the Cauchy problem for the one-dimensional heat equation. We consider
ut = uxx ,
x (, ) ,
t (0, T ] ,
(5.7)
(5.8)
This can be treated in a manner analogous to what was done for ODE initial-value problems in
Chap. 4.
Recall that for such problems, the easiest scheme to construct is the forward-Euler method,
which consists of replacing the time derivative with a forward difference, and evaluating the righthand side at the previous time level. Thus, (5.7) becomes
n
un+1
m um
= (uxx )nm
k
or
n
n
un+1
m = um + k(uxx )m .
(5.9)
In (5.9) the m-subscripts denote spatial grid points, while the n-superscripts indicate the time level.
To complete the discretization we replace uxx with a centered-difference approximation and obtain
n
un+1
m = um +
k
n
n
n
u
2u
+
u
,
m1
m
m+1
h2
m = 1, 2, . . . , M .
(5.10)
As indicated in Fig. 5.2, k denotes the time step, and h is the space step. Recall that in the context
of ODEs, Eulers method is a single-step procedure; in the PDE context we use the terminology
two-level scheme since contributions from two time levels, n and n + 1, appear in the discrete
formula.
It is easily seen that (5.10) is an explicit method, since the new time level values at each grid
point can be directly calculated from previous ones. However, it should be apparent that we cannot
compute all values at time level n + 1 from (5.10) using only the values given at time level n because
at each of the endpoints of the computational interval the right-hand side of (5.10) requires one
additional point. This is indicated in Fig. 5.2. Here, we show the so-called grid stencil, or mesh
star, for the general mth and final M th grid points. The dashed line in the latter stencil indicates
that there is no grid point value available for evaluation of that part of the stencil.
There are several methods by which this situation can be remedied; we will describe only one
of them here. It is to calculate only between m = 2 and M 1 at the second time level, between
125
n+1
k
n
...
m1
m+1
...
M1
Figure 5.2: Mesh star for forward-Euler method applied to heat equation
m = 3 and M 2 at the third, etc. That is, we lose one spatial grid point from each end of the
discrete domain with each new calculated time step. Thus, we must start with a sufficient number
of grid function values at the initial time to guarantee that we cover the desired spatial interval at
the final time. This approach is computationally inefficient because more grid points are computed
than are needed in all but the final step. On the other hand, stability analyses (see below) for such
a method are much more straightforward than for alternatives employing so-called outflow or
nonreflecting boundary conditions.
We now provide a pseudo-language algorithm for implementing this method.
Algorithm 5.1 (Forward-Euler/Centered-Difference Approximation for 1-D Heat Equation)
1. Input number of time steps Nt , number of space steps at final time Nx , time step size k, and
[ax , bx ] the solution domain at final time.
2. Input initial data: u0i , i = 1, 2, . . . , Nx .
3. Load zeros at additional points needed to obtain solution on required final interval:
Do i = 1, Nt
u01i = 0
u0Nx +i = 0
Repeat i
4. Calculate spatial step size, h: h = (bx ax )/(Nx 1).
5. Begin time stepping:
Do n = 1, Nt
mstrt = n Nt
mstop = Nx + Nt n
Do m = mstrt, mstop
n1 + k un1 2un1 + un1
unm = um
m
m+1
m1
h2
Repeat m
Output results for time step n
Repeat n
126
The main reason for introducing the explicit approximation (5.10) is that it is quite simple,
and it provides an easy starting point for the discussions of consistency and stability of difference
approximations to partial differential equations. By consistency we mean that the difference approximation converges to the PDE as h, k 0, and by stability we mean that the solution to the
difference equation does not increase with time at a faster rate than does the solution to the differential equation. Consistency and stability are crucial properties for a difference scheme because
of the following theorem due to Lax (see Richtmyer and Morton [24]).
Theorem 5.1 (Lax Equivalence Theorem) Given a well-posed linear initial-value problem, and
a corresponding consistent difference approximation, the resulting grid functions converge to the
solution of the differential equation(s) as h, k 0 if and only if the difference approximation is
stable.
It is important to understand the content of this theorem. One can first deduce the less than
obvious fact that consistency of a difference approximation is not a sufficient condition for guaranteeing that the grid functions produced by the scheme actually converge to the solution of the
original differential equation as discretization step sizes are refined. In particular, both consistency
and stability are required. As will be evident in what follows, consistency of a difference approximation is usually very straightforward though sometimes rather tedious to prove, while proof of
stability can often be quite difficult.
Consistency of Forward-Euler/Centered-Difference Approximation to
the Heat Equation
We will now demonstrate consistency of the difference scheme (5.10). As a byproduct we also
obtain the order of accuracy of the difference approximation. For this analysis it is convenient to
rewrite (5.10) as
k
2k
n+1
(5.11)
um 1 2 unm 2 unm1 + unm+1 = 0 .
h
h
k2
(utt )nm +
2
h3
h2
= unm h(ux )nm + (uxx )nm (uxxx )nm +
2
6
2
3
h
h
= unm + h(ux )nm + (uxx )nm + (uxxx )nm +
2
6
un+1
= unm + k(ut )nm +
m
unm1
unm+1
h4
(uxxxx )nm
24
h4
(uxxxx )nm + ,
24
and if we substitute these expansions into (5.11) and rearrange the result, we find
(ut
uxx )nm
k
h2
n
n
+
(utt )m (uxxxx )m + = 0 .
2
12
(5.12)
The first term in this expression is the original differential equation evaluated at the arbitrary
space-time point (xm , tn ). Thus, the difference approximation converges to the differential equation
provided the second term 0 with h, k (under the assumption that all neglected terms are indeed
small compared with those retained). This will usually hold (but somewhat pathological exceptions
can be constructed) whenever utt and uxxxx are continuous functions of (x, t); i.e., u C 2,4 . Hence,
consistency has been shown within this context.
127
The second term of (5.12) is the dominant, or principal, term of the global truncation error.
This is analogous to the ODE case discussed
in Chap. 4. In order notation we also have that the
local truncation error isO k k + h2 . The global (in time) truncation error, and thus the order of
the method is O k + h2 . Hence, we have shown that the difference scheme based on forward-Euler
time integration and centered spatial differencing is consistent with the heat equation, and it is
first order in time and second order in space. This is the expected result.
Stability of Forward-Euler Approximation to the Heat Equation
We next carry out a stability analysis for the difference equation (5.10). Define r k/h2 , and now
write (5.10) as
un+1
= (1 2r)unm + r unm1 + unm+1 .
(5.13)
m
It is not altogether obvious that difference equations can be unstable, particularly when they are
consistent approximations to a well-posed PDE problem, but recall the situation for ODE IVPs
discussed in Chap. 4; we will here demonstrate that analogous behavior occurs in the PDE case.
It is worth noting that were this not true the consequences of the Lax theorem would be mere
triviality.
We have already seen that for a consistent difference approximation the difference between the
theoretical solution of the PDE and the theoretical solution of the difference equation is related
to the truncation error, and hence for sufficiently small h, k is always bounded. However, the
theoretical solution to the difference equation will not be obtained in actual calculations, mainly
because of machine rounding errors; so we must define a new error component as the difference
n of the difference equation: we
between the theoretical solution unm and the computed solution vm
set
n
n
zm
= unm vm
,
n is the result of actual computation, again analogous to the ODE IVP treatment of Chap.
where vm
n as n increases, for if this quantity grows rapidly it
4. We wish to investigate the growth of zm
follows that the computed solution may bear no resemblance to the true solution, and in such a
case Eq. (5.13) will be called unstable.
To proceed with this investigation we note, as we have already seen in the ODE case, that for
n must satisfy the same difference equation as does un and v n . Hence,
linear difference equations zm
m
m
in the present case, from (5.13) we have
n+1
n
n
n
zm
= (1 2r)zm
+ r(zm1
+ zm+1
).
(5.14)
n }, has a Fourier
We will now assume that z(x, t), the continuous analog of the grid function {zm
representation
X
z(x, t) =
Al (t)eil x
l
for arbitrary wavenumbers l . We observe that this is a rather natural assumption in light of the
form of exact solutions to the heat equation (see, e.g., [2]). It is sufficient to consider a single
arbitrary term of this series because if any one term is growing unboundedly, then this must also
be true, in general, for the series.
Now we want the initial error at t = 0 to consist only of error in the spatial representation, so
we set
n
(5.15)
zm
= et eixm = enk eimh = n eimh ,
128
where the first right-hand equality follows from exact solutions to the heat equation but now with
0 = eimh ; furthermore, this will not grow
the l index suppressed. Clearly, at t = 0, the error is zm
in time, provided
|| 1 .
(5.16)
This is the well-known von Neumann stability criterion for stability of discretizations of parabolic
PDEs.
For two-level methods, and only a single differential equation, it is both necessary and sufficient
for stability; but for schemes employing three (or more) time levels, and/or applied to more than
one equation, it is only necessary, but not sufficient in general for stability.
We now substitute the right-hand side of (5.15) into (5.14) to obtain
n+1 eimh = (1 2r) n eimh + r n ei(m1)h + n ei(m+1)h .
Division of n eimh leaves
= (1 2r) + r eih + eih
= 1 2r(1 cos h)
h
= 1 4r sin2
.
2
Thus, in order for the forward-Euler method to satisfy the von Neumann condition (5.16) we must
choose r so that
h
1
h .
1 1 4r sin2
2
It is easily seen that the right-hand inequality is satisfied r 0, while the left-hand side implies
r
1
2 sin2
h
2
We must choose h so that the right-hand side is a minimum in order to place the strongest restrictions on r. This minimum occurs when the denominator is a maximum, namely when sin2 h
2 = 1.
Thus, the stability requirement for the forward-Euler scheme is
0r
1
.
2
(5.17)
(Of course we would never take r = 0 because this implies a zero time step.)
If we recall that r = k/h2 , we see an important consequence of this stability criterion: it is
that as the spatial grid is refined to improve accuracy, the time steps must be decreased as the
square of the spatial grid size to maintain stability. For example, if h = 0.1, then we must have
k 0.005. This is often a far too restrictive requirement, and it has motivated development of
implicit methods such as those that we next consider.
5.3.2
Although the backward-Euler method is only first-order accurate in time, it is nevertheless widely
used because of its favorable stability properties and its simplicity. Recall from our discussions
for ODEs that the backward-Euler procedure is A-stable. We will see that for the PDE case it is
unconditionally stable for linear parabolic problems.
129
We now consider the heat equation on a bounded domain; i.e., we treat the initial boundary
value problem
ut = uxx , x (0, 1) , t (0, T ] ,
(5.18)
with initial conditions
u(x, 0) = u0 (x) ,
(5.19)
u(1, t) = g(t) .
(5.20)
Here u0 (x), f (t) and g(t) are known functions. To approximate (5.18) with the backward-Euler
method, we employ a backward difference for the time derivative approximation, and then in
contrast to (5.10), evaluate everything on the right-hand side at the advanced time level. Hence,
n1
unm um
= (uxx )nm ,
k
or after shifting temporal indexing forward by one and introducing a centered-difference spatial
approximation for uxx ,
k n+1
n+1
n
n+1
un+1
(5.21)
m = um + 2 (um1 2um + um+1 ) .
h
Now unknown n + 1 time level results appear on both sides of the equation, so it is not possible to
explicitly solve for un+1
m ; hence, this is an implicit, but still two-level, method.
As is standard in treating such difference equations, we first move all n + 1 time level terms to
the left-hand side and rearrange this as
1
1 n
n+1
n+1
um1 +
+ 2 un+1
m = 2, 3, . . . , M 1 .
(5.22)
m um+1 = um ,
r
r
The collection of all such equations results in a tridiagonal system for the grid function {un+1
m }.
Clearly, for m = 1 we have
un+1
= f n+1 ,
1
and for m = M
n+1
.
un+1
M =g
Just as was the case for ODE BVPs, these boundary conditions can be inserted directly into the
n+1
and un+1
system for un+1
m , or they can be used to eliminate u1
M from the equations corresponding,
respectively, to m = 2 and m = M 1 in (5.22). (Recall Fig. 4.9.)
We will not present a detailed pseudo-language algorithm for the backward-Euler method at
this time because it is very similar to that for the trapezoidal method to be presented later and, in
fact, can be embedded as an option in this latter algorithm.
Consistency and Stability of Backward-Euler Method
Consistency (and order of accuracy) of the backward-Euler method can be demonstrated in precisely
the same way as already carried out for the forward-Euler discretization, and the result is the same.
Namely, as we should expect, backward Euler/centered-discretizations are (globally) first-order
accurate in time and second order in space. We leave demonstration of this as an exercise for the
reader.
The von Neumann stability analysis applied to (5.22) shows that the method is stable for any
value of r > 0; i.e., it is unconditionally stable. We leave this also as an exercise to the reader.
130
It should be noted, however, that such an analysis does not include the effects of either boundary
conditions or variable coefficients. Thus, we very briefly consider a method which does, but one
that is rather seldom used because it, in general, requires use of additional numerical methods to
obtain sharp predictions regarding stability.
The method is known as the matrix method, or matrix stability analysis, and it is based on the
fact that any two-level method for PDEs can be formally represented as
An+1 un+1 = B n un + f n ,
(5.23)
where A and B are both square matrices whose elements depend upon the specific difference
approximation employed (and, of course, on the specific PDE being solved), and the superscripts
indicate that time-dependent coefficients can be treated; f n is a given nonhomogeneous function
and/or boundary data. (Recall the manner in which boundary conditions entered the matrix
structure of finite-difference approximations to ODE BVPs in Chap. 4.) It is of interest to observe
that if the scheme is explicit, An+1 = I n. By the same arguments employed in the von Neumann
stability analysis, we can show for the linear case that the error vector z n satisfies (5.23) without
the inhomogeneous term. We have
An+1 z n+1 = B n z n ,
or after formally solving for z n+1 ,
z n+1 = An+1
1
B nzn .
(5.24)
where k k denotes the spectral norm. Checking this requires computation of the largest eigenvalue
h
1 n i h n+1 1 n iT
of An+1
A
B
, which can seldom be done analytically, and is usually rather
B
expensive if done numerically. Nevertheless, this approach does have some theoretical applications,
as we shall see later, and with ever-increasing hardware performance it is beginning to represent a
tractable approach for stability analysis for a wide range of numerical PDE problems.
5.3.3
To this point, the (two) methods considered have been only first-order accurate in time. Often,
we wish to have more accuracy to permit taking larger time steps, thus reducing rounding errors,
or to obtain better accuracy for any given time step size. Hence, we consider two second-order
methods. The first is a famous method, due to Richardson, which is consistent with the heat
equation to second order in both space and time, but which is unconditionally unstable, and thus
useless as a solution method. The second is probably the most widely-used scheme for solving
parabolic equations. It is based on trapezoidal integration; it is also second order in both space
and time, and is unconditionally stable for linear constant coefficient problems. It is known as the
CrankNicolson method.
Richardsons Method for the Heat Equation
Richardsons scheme is the most natural second-order method. It is constructed by replacing all
derivatives, in both space and time, with second-order centered differences. Thus, for the heat
131
or
n1
un+1
m = um +
2k n
(u
2unm + unm+1 ) .
h2 m1
(5.25)
We observe that this is a three-level difference scheme, and as a consequence it can have more than
one solution. Moreover, because it is centered in time it has a backward as well as a forward (in
time) solution, and this is the underlying reason for its instability. Recall that the backward heat
equation problem is ill posed, and in general has no bounded solution. This lack of boundedness in
the true solution leads to lack of boundedness in the solution to the difference equation as well for
any consistent discretization admitting multiple solutions which can propagate backwards in time.
We note that a von Neumann stability analysis can be applied to (5.25); however, it is not as direct
as in the case of two-level schemes (it depends on conversion of the multi-level method to a system
of two-level difference equations; see e.g., Richtmyer and Morton [24] and treatment of the wave
equation in Sec. 5.5 below). Furthermore, the results obtained provide necessary, but not generally
sufficient, stability criteria.
We must again emphasize that this second-order discretization should never be used for solving
the heat equation, or any other parabolic equation. On the other hand, if one is merely seeking
a numerical evaluation of the heat equation (for example, to confirm heat balances in laboratory
experiments) at a fixed instant in time, then Eq. (5.25) is a quite satisfactory approximation. The
difference between this type of application and using (5.25) to solve the heat equation is that
evaluation does not imply evolution, and it is (time-like) evolution of difference approximations
that can lead to instability.
CrankNicolson Method for the Heat Equation
To derive the CrankNicolson method we will employ a procedure that has received wide acceptance
in the study of difference approximations of evolution, i.e., time-dependent, equations. Namely, we
view the PDE as an ODE in time. Thus, we can write the heat equation as
du
= uxx F (u) ,
dt
(5.26)
with initial and boundary conditions given, for example, by Eqs. (5.19), (5.20) and develop solution
procedures based on ODE methods. It is worth mentioning at this point that it can be shown that
the eigenvalues of F (u) are negative and large in magnitude. Hence (5.26), viewed as a system of
ODEs, is stiff, and we should seek an A-stable method for its solution.
Thus, we apply the trapezoidal rule to (5.26) and obtain
un+1 = un +
k
k
F (un+1 ) + F (un ) = un + (uxx )n+1 + (uxx )n .
2
2
We already know that this is (globally) second order in time, and it will be second order in space
if we approximate uxx with a second-order centered difference. Hence, we write
un+1
= unm +
m
k n+1
n+1
n
n
n
(um1 2un+1
m + um+1 ) + (um1 2um + um+1 ) ,
2
2h
132
m = 2, 3, . . . , M 1 . (5.27)
Since the right-hand side consists only of already known grid-function values at time level n, this
constitutes a tridiagonal system for the implicit determination of the solution vector un+1 at time
level n + 1.
We should comment that division by r to obtain Eq. (5.27) is rather customary (although the
factor of 21 is usually not included in the definition of r) and analogous to multiplying by h2 in the
ODE case. Also recall, as observed in Chap. 1, that pivoting is seldom performed in the context of
tridiagonal LU decomposition, and that diagonal dominance is required to control round-off errors.
Clearly, as r 0 the left-hand side of (5.27) will become strongly diagonally dominant.
It is of interest to analyze the consistency and stability of this widely-used method. In the
process we will obtain the formal truncation error. We note at the outset that (5.27) has been
arranged in a different manner than was the case for the explicit Euler formula (5.11). In particular,
to obtain (5.27) we have multiplied by 1/r. Our results will thus be global in time. We also carry
out the analysis in terms of the (n + 1)th time level, rather than the nth as done previously; but
this is an unimportant detail that cannot effect the results. (We leave as an exercise for the reader
to repeat the following analysis for time level n.) To begin we write (5.27) as
1
1
n+1
n+1
n
n
um 2
unm + un+1
(5.28)
um1 + um1 2 +
m+1 + um+1 = 0 .
r
r
We expand unm1 about the n + 1th time level to obtain
k3
k2
(utt )n+1
(uttt )n+1
m1
m1 + .
2
6
Next, we expand each of the terms appearing on the right in terms of the mth spatial grid point.
Thus,
n+1
n+1
n
un+1
m1 + um1 = 2um1 k(ut )m1 +
h3
h4
h2
(uxx )n+1
(uxxx )n+1
(uxxxx )n+1
m
m +
m
2
6
24
h2
h3
n+1
= (ut )n+1
(utxx )n+1
(utxxx )n+1
m h(utx )m +
m
m +
2
6
h2
n+1
(uttxx )n+1
= (utt )n+1
m
m h(uttx )m +
2
n+1
= (uttt )n+1
m h(utttx )m +
n+1
n+1
un+1
m1 = um h(ux )m +
(ut )n+1
m1
(utt )n+1
m1
(uttt )n+1
m1
133
6
In addition, we have
1
k3
1
1
k2
n
2+
u 2
um = 4u + 2
kut utt + uttt .
r
r
r
2
6
We now combine these results and use the definition of r to obtain, after some rearrangement,
ut uxx
k2
k2
k
h2
uxxxx + uttt uttxx + (utxx utt ) + = 0 .
12
6
4
2
Finally, we observe that since ut = uxx , it follows that utt = uxxt , and if u is sufficiently smooth
uxxt = utxx . Hence, in this case the O(k) term is identically zero, and since uttxx = uttt , the O(k2 )
terms can be combined so that the above becomes
(ut uxx )n+1
=
m
1
2
n+1 2
(uxxxx )n+1
O(h2 + k2 ) .
m h + (uttt )m k
12
(5.29)
The right-hand side of Eq. (5.29) is the global truncation error; it clearly goes to zero as h, k 0,
proving consistency of the CrankNicolson scheme. It should be clear at this point that our ODE
analogy breaks down to some extent because of the additive nature of the different portions of
the truncation error. In particular, letting only k 0, but with h finite, will not result in the
exact solution to a PDE problem unless the solution happens to be such that uxxxx and all higher
derivatives are identically zero (as, for example, with a low-degree polynomial for the true solution).
In fact, we would obtain the exact solution to a different problem,
ut = uxx + T (h) ,
where T (h) is the spatial truncation error due to a step size h. It is important to keep this point in
mind when performing grid-function convergence tests for PDEs because if h remains fixed while
k is decreased, convergence will occur, but not to the true solution of the problem. This can often
lead to seemingly ambiguous results. For a method with truncation error O(h2 +k2 ), k and h should
be simultaneously reduced by the same factor during convergence tests. Analogous arguments and
conclusions apply if h 0 with k finite.
We next briefly consider a stability analysis for the CrankNicolson procedure. We would
expect, because trapezoidal integration is A-stable, that the CrankNicolson scheme should be
unconditionally stable for linear constant-coefficient problems. This is, in fact, the case as we will
show. We apply the von Neumann analysis discussed earlier, but once again note that this does
not account for boundary conditions. Again letting
n
= n eimh
zm
134
denote the difference between the true solution and the computed solution to (5.27), we obtain
=
2 1r 2 cos h
.
2 + 1r + 2 cos h
Recalling that the von Neumann stability condition (which is both necessary and sufficient in this
case) is || 1, we see that the following two inequalities must be satisfied:
2 1r 2 cos h
1,
2 + 1r + 2 cos h
and
2 1r 2 cos h
1.
2 + 1r 2 cos h
It is a simple exercise, left to the reader, to show that both of these hold, completely independent of the value of r. Hence, the CrankNicolson method applied to the linear heat equation is
unconditionally stable.
Finally, we remark that this result is not significantly altered by the presence of inhomogeneities;
however, in such cases we employ a more general concept of stability. Namely, because the true
solutions may grow as fast as exponentially with time, it is too restrictive to require || 1. Instead,
generally for nonhomogeneous problems, we only require
|| 1 + O(k)
for stability.
We conclude this section on second-order methods for parabolic equations with a pseudolanguage algorithm for the trapezoidal method. Before presenting this we provide a slightly generalized version of Eq. (5.27) allowing us to write a single code that embodies all of forward- and
backward-Euler, and trapezoidal, methods for initial boundary value problems for the heat equation. This is done by viewing the construction of the trapezoidal method as consisting of a simple
average of values between the n and n + 1 time levels, and replacing this with a weighted average. If
we denote this weight by , (0 1), then the basic trapezoidal integration formula is replaced
by the more general form
k
n+1
n
n+1
n
+ (1 ) (um1 2um + um+1 )
.
um = um + 2 (um1 2um + um+1 )
h
It is easily seen that = 21 yields the usual trapezoidal formula. If we set = 1, we obtain the
backward-Euler method, while = 0 results in the forward-Euler approximation.
One often encounters the terminology semi-implicit, implicit and explicit for these three
cases, but the first of these is misleading because as will be evident from the pseudo-languange algorithm that follows, the trapezoidal method is no less implicit than is the backward-Euler method:
the entire spatial grid function is computed simultaneously (implicitly) at each time step in both
approaches. They do, however, exhibit different behavior for large values of time step k. In particular, despite the fact that trapezoidal integration is A-stable, if k is too large solutions will
oscillate (boundedly) from one time step to the next. For this reason backward Euler is sometimes
preferred even though it is only first-order accurate in time. An instance of this is the case of
calculating a steady-state solution by integrating a formally time-dependent equation to steady
state, a widely-used practice in computational fluid dynamics and heat transfer.
If we define
(1 )k
k
r2
,
r1 2 ,
h
h2
135
um + um+1 = um1 +
unm unm+1 .
um1 2 +
r1
r1
r1
r1
r1
(5.30)
We now present the pseudo-language algorithm for implementing this generalized trapezoidal
method.
Algorithm 5.2 (Generalized Trapezoidal Method)
1. Input number of time steps Nt , number of spatial grid points Nx , time step size k, and
endpoints of spatial domain, ax , bx , time integration weight, .
2. Input initial data, u0m , m = 1, . . . , Nx , and initial time, t0 .
3. Calculate grid spacing: h = (bx ax )/(Nx 1), time integration parameters, r1 = k/h2 , r2 =
(1 ) k/h2 .
4. Begin time stepping
Do n = 1, Nt
tn = k n + t0
2r2
n1
n1
+ um+1
+ r1
Bm = rr12 um1
Repeat m
A1,Nx = 0
A2,Nx = 1
A3,Nx = 0
BNx = Rhtbndry(tn )
1
r1
n1
um
136
5.3.4
Our final topic in this study of parabolic equations is treatment of the two-dimensional heat equation. We first point out that, generally, explicit methods should not be used because the stability
requirements are even more stringent for 2-D algorithms than in the 1-D case, leading to usually
unacceptable increases in required arithmetic to integrate to a desired final time. Thus, we will
restrict attention to implicit methods. In fact, we shall consider only the 2-D analogue of the
CrankNicolson method although a similar approach based on backward-Euler time integration is
also widely used.
The problem we study is the 2-D heat equation on the unit square,
ut = uxx + uyy + f (x, y, t) ,
t (0, T ] ,
(5.31)
with prescribed Dirichlet conditions on the boundaries, , and given initial data on .
If we apply the trapezoidal integration to (5.31) we obtain for interior grid points
n
un+1
l,m = ul,m +
i
kh
n+1
3
n
n
+
f
+
f
(uxx + uyy )n+1
+
(u
+
u
)
xx
yy l,m
l,m + O(k ) ,
l,m
l,m
2
(5.32)
un+1
,
+
f
f
+
=
u
+
(u
+
u
)
(u
+
u
)
xx
yy l,m
xx
yy l,m
l,m
l,m
l,m
2
2
2 l,m
and then replace spatial derivatives with centered differences, constructed with uniform grid spacing
h in both x and y directions, to obtain
un+1
l,m
k n+1
n+1
n+1
n+1
n+1
u
+
u
4u
+
u
+
u
l1,m
l,m1
l,m
l,m+1
l+1,m
2h2
k n+1
k
n
). (5.33)
+ fl,m
= unl,m + 2 unl1,m + unl,m1 4unl,m + unl,m+1 + unl+1,m + (fl,m
2h
2
If we suppose that the spatial grid consists of Nx points in the x direction, and Ny in the y
direction (with Nx = Ny so that hx = hy = h), then (5.33) holds for all l = 2, 3, . . . , Nx 1 and
m = 2, 3, . . . , Ny 1. For a Dirichlet problem, if l = 1 or Nx , or m = 1 or Ny , we replace (5.33)
with the appropriate Dirichlet condition,
ul,m = (xl , ym , t) ,
(5.34)
137
where (x, y, t) is a prescribed boundary function. We note that boundary condition evaluation for
the left-hand side of (5.33) must be at time level n + 1, while that on the right-hand side must be at
time level n. Failure to adhere to this will result in degradation of formaland observedaccuracy.
Before writing the final form of the difference equations we again introduce the notation r
k/2h2 . We now divide by r in (5.33) and obtain
un+1
l1,m
un+1
l,m1
1
4+
r
n+1
n+1
un+1
l,m + ul,m+1 + ul+1,m
1
n
n
n
n+1
= ul1,m ul,m1 + 4
).
+ fl,m
unl,m unl,m+1 unl+1,m h2 (fl,m
r
We observe that everything on the right-hand side is known prior to the start of calculations for
time level n + 1, so we define
1
n+1
n
n
n
n
Fl,m ul1,m ul,m1 + 4
).
+ fl,m
unl,m unl,m+1 unl+1,m h2 (fl,m
r
Hence, the difference equations at interior grid points can be expressed simply as
1
n+1
n+1
n+1
n+1
n
ul1,m + ul,m1 4 +
un+1
l,m + ul,m+1 + ul+1,m = Fl,m ,
r
(5.35)
(5.36)
n+1
Moreover, (5.34) is also of this form if we set A3l,m = 1, Ail,m = 0, i 6= 3, and bl,m = l,m
. The
complete matrix associated with this system of difference equations has the structure displayed in
Fig. 5.3.
This matrix consists of only five nonzero bands. However, the matrix structure shown is such
that this high degree of sparsity cannot be readily exploited. Thus, usual direct elimination cannot
be applied. It is the presence of the two outer bands, symmetrically spaced Ny elements from the
main diagonal, that prevents use of sparse LU decomposition. It turns out, however, that with a
little manipulation, the original system of difference equations can be converted to a system of two
sets of equations, each of which involves only one spatial direction, and requires solution of only
tridiagonal systems. Such procedures are described generically as alternating-direction-implicit
(ADI) methods. They require only O(N ) (N = Nx Ny ) arithmetic operations per time step, in
contrast to the O N 3 required by direct Gaussian elimination applied to the entire system.
To derive the ADI scheme to be used here, we observe that the approximation (5.33) can be
rearranged and represented as
k 2
k n+1
k 2
n+1
2
2
n
+ fl,m
)
(5.37)
I (D0,x + D0,y ) ul,m = I + (D0,x + D0,y ) unl,m + (fl,m
2
2
2
2
2 u
l = 2, . . . , Nx 1, m = 2, . . . , Ny 1. Now each of D0,x
l,m and D0,y ul,m is just a tridiagonal
matrix times the vector grid function {ul,m }, as we have already seen in the study of ODE BVPs.
For notational simplicity, label these matrices as Ax and Ay , respectively, and write the above as
k
k
k
n+1
= I + (Ax + Ay ) un + (f n+1 + f n ) .
(5.38)
I (Ax + Ay ) u
2
2
2
138
Ny + 1
b1,1
u1,1
u1,2
.
.
.
u1, Ny
b1,2
.
.
.
b 1, Ny
u2,1
.
.
.
ul,m
.
.
.
.
uNx ,1
.
.
.
.
uNx ,Ny
b2,1
.
.
.
b l,m
.
.
.
.
bNx ,1
.
.
.
.
bNx ,Ny
l, m
139
the CrankNicolson procedure with a system involving only tridiagonal matrices has not yet been
achieved because (5.39) contains matrix productsnot single tridiagonal matrices.
The standard remedy is to split (5.39) into two equations as follows:
k
k n
k
n+1
n
A
A
(5.40a)
u
=
I
+
I
x
y u + f ,
2
2
2
k
k
k
n+1
I Ay u
= I + Ax un+1 + f n+1 .
(5.40b)
2
2
2
Thus, we calculate the advanced time step grid-function values in two consecutive steps. The first
involves only x derivatives at the advanced time level (i.e., at n + 1 ), while the second involves
only y derivatives. In each case the advanced time values are calculated implicitly, thus giving rise
to the name alternating direction implicit, or ADI. There are many such procedures, and the one
presented here is among the oldest, but still widely used, due to Peaceman and Rachford [23].
1
It is rather common practice to view un+1 as being equivalent to un+ 2 . In general, this is
not true for ADI-like schemes. In particular, for problems involving time-dependent boundary
conditions, second-order temporal accuracy can be lost if boundary conditions for the first split
1
step are evaluated at t = tn+ 2 . The correct implementation of time-dependent Dirichlet conditions
is presented in Mitchell and Griffiths [21]; we shall not pursue this topic further herein except
to mention that in the specific case of PeacemanRachford ADI presented here, numerical tests
1
k
I Ax un+1 = I + Ay un ,
(5.41a)
2
2
k
k
n+1
(5.41b)
= I + Ax un+1 .
I Ay u
2
2
If we solve the first of these for un+1 , and substitute the result into the second, we obtain
1
k
k
k
k
n+1
I Ay u
= I + Ax
I Ax
I + Ay un .
2
2
2
2
We now multiply by I k2 Ax , and note that the matrices I + k2 Ax and I k2 Ax commute,
leading to
k
k
k
k
n+1
I Ay u
= I + Ax
I + Ay un .
I Ax
2
2
2
2
This is exactly the unsplit scheme, (5.39), with f 0. Hence, for constant coefficient homogeneous
problems, second-order global accuracy is achieved up to implementation of boundary conditions.
In fact, this is still true for f 6 0 (and even time-dependent); but this is somewhat more tedious
to demonstrate, and we leave this to the ambitious reader.
Stability of PeacemanRachford ADI
We next consider stability of the PeacemanRachford scheme. We shall not carry out a rigorous
analysis, but instead give an heuristic treatment that is somewhat specific to the 2-D heat equation.
140
For this particular case, it is easily seen that Ax = Ay up to multiplication by a permutation matrix.
We then recognize that each step of (5.41) corresponds to a 1-D CrankNicolson procedure, which
previously was shown to be unconditionally stable in the sense of von Neumann. Hence, if we
n defined as in (5.15),
ignore boundary conditions, we can infer that for zl,m
n
n+1
,
zl,m = |x | zl,m
and we have |x | 1 for the x-direction step. Similarly, for the y direction we have
n+1
n+1
zl,m = |y | zl,m
with |y | 1 . Thus, it follows that
n
n+1
zl,m = |y | |x | zl,m
n
,
= || zl,m
k
2h2
(5.42)
(5.43)
We remark that a similar treatment can be carried out when the grid spacing h is different in the
two different directions, say hx and hy , but the resulting equations are considerably more tedious.
We leave analysis of this as an exercise for the reader.
141
There are several things to notice regarding Eqs. (5.42) and (5.43). The first is that each of
these separate steps is formally very similar to the 1-D CrankNicolson scheme, a fact we have already used to make plausible the unconditional stability of the PeacemanRachford procedure. The
second is that (5.42) holds for each m = 1, 2, . . . , Ny except when Dirichlet conditions are imposed
at one, or both, y-direction boundaries. In particular, for each value of
index m, (5.42) results
the
n+1
n+1
n+1
in a tridiagonal system of Nx equations for the grid function values u1,m , u2,m
,...,
, . . . , ul,m
o
n
n+1
n+1
we use O (Nx Ny ) O(N ) arithuNx ,m . Hence, to compute the complete grid function ul,m
metic operations.
Construction of the right-hand sides of (5.42) deserves some special attention. We first note
that all quantities are known either from previous computation or from given data; but some care
must be exercised in loading results from previous calculations. In particular, we must be sure to
use only results from time level n and not recently-calculated grid-function values from n + 1 in
constructing the right-hand side of (5.42). Figure 5.4 indicates why this is a concern. At the mth
Boundary #4
Boundary #3
m+1
m
m 1
. . . .
Boundary #1
. . . .
Ny
2
1
1
. . . .
. . . .
Nx
Boundary #2
142
conclude that the total arithmetic per time step is O (Nx Ny ), and since N = Nx Ny is the length of
the solution vector, the PeacemanRachford scheme requires only O(N ) arithmetic operations per
time step.
We now present a pseudo-language algorithm for implementing the PeacemanRachford scheme.
Algorithm 5.3 (PeacemanRachford ADI Method for Dirichlet Problems)
1. Read inputs: Nt , Nx , Ny , k, h, and initial data
2. Calculate r =
k
2h2 ;
3. Print inputs
4. Increment time: tn = t0 + n k
5. Calculate n + 1 intermediate time level results
Do m = 1, Ny
If m = 1 then
Do l = 1, Nx
usv,l = unl,m
ul,m = Bndry xl , ym , tn + k2 , 2
Repeat l
Else if m = Ny
Do l = 1, Nx
ul,m = Bndry xl , ym , tn + k2 , 4
Repeat l
Else
Do l = 1, Nx
If l = 1 or l = Nx then
Else
A1,l = 0
A2 , l = 1
A3 , l = 0
If l = 1, B1 = Bndry xl , ym , tn + k2 , 1
If l = Nx , BNx = Bndry xl , ym , tn + k2 , 3
A1,l = 1
A2,l = 2 + 1r
A3,l = 1
Bl = usv,l + 2 1r unl,m unl,m+1 h2 f n (xl , ym , tn )
End if
Repeat l
Call LUDCMP (A, B, Nx )
Do l = 1, Nx
usv,l = unl,m
n+1
= Bl
ul,m
End if
Repeat m
Call LUDCMP (A, B, Ny )
Do m = 1, Ny
n+1
usv,m = ul,m
un+1
l,m = Bm
Repeat m
End if
Repeat l
Print results for t = tn+1
7. Increment n
143
144
Although the preceding algorithm has been constructed specifically for the heat equation, the
extension to more general 2-D parabolic problems is straightforward. On the other hand, the
PeacemanRachford procedure cannot be implemented for 3-D problems in any direct manner.
Nevertheless, there are numerous related methods that can be applied in this case. We refer the
reader to Douglas and Gunn [7] for a very general treatment.
In closing this section on the PeacemanRachford ADI method we should point out that in
recent years iterative techniques have often been employed to solve the five-band system shown in
Fig. 5.3, especially in FEMs. It is important to recognize, however, that such approaches cannot
possibly compete with the PeacemanRachford scheme (or other similar time-splitting methods)
in terms of total required arithmetic. Moreover, ADI schemes in general are highly parallelizable,
so they present several distinct advantages on modern multiprocessors that are difficult to realize
with various alternative approaches.
5.4
Elliptic Equations
Recall that, as discussed earlier, constant coefficient elliptic PDEs in two space dimensions can
always be reduced to the canonical form
uxx + uyy + Au = f (x, y) ,
(5.44)
where A = 1, 0 or +1. In the present treatment we will consider only the case A = 0, so that
(5.44) becomes a Poisson equation. Moreover, we will restrict attention to Dirichlet problems on
rectangular domains. We remark that the cases A = 1 can be treated completely analogously to
that considered here, and other types of boundary conditions can be implemented as described in
Chap. 4 for ODE problems.
Two classical methods will be applied to the solution of
uxx + uyy = f (x, y) ,
(5.45)
with
u(x, y) = g(x, y) ,
(x, y) .
(5.46)
These methods are: i) successive overrelaxation (SOR) and ii) alternating direction implicit (ADI).
We have already encountered forms of both of these methods in our earlier studies. It must be
noted that although many other more efficient procedures for solving elliptic problems are now in
existence, SOR and ADI still are widely used, and understanding them helps to form the basis for
developing more sophisticated modern methods.
We shall begin by discretizing (5.44) in our usual manner, viz., by replacing derivatives with
second-order centered differences. Thus, at grid point (i, j) we have
ui,j1 2ui,j + ui,j+1
ui1,j 2ui,j + ui+1,j
+
= fij ,
2
hx
h2y
or (with hx = hy = h)
ui1,j + ui,j1 4ui,j + ui+1,j + ui,j+1 = h2 fij
(5.47)
145
i = 1, Nx with j = 1, 2, . . . , Ny ,
j = 1, Ny with i = 1, 2, . . . , Nx .
(5.48)
This set of equations can be formally solved with the system (5.47), or it can be used to eliminate
terms from those equations in (5.47) corresponding to points adjacent to boundaries. For general
codes, written to permit other boundary condition types as well as Dirichlet conditions, the former
is preferred, as has been noted previously.
We leave formal truncation error analysis of this discretization as an exercise for the reader.
Our intuition should suggest (correctly) that the dominant truncation error is O(h2 ) for sufficientlysmooth solutions u(x, y). Moreover, because the problem being considered is independent of time,
and correspondingly the difference approximation is not of evolution type, stability of the approximation is not an issue.
The matrix form of (5.47), (5.48) is presented in Fig. 5.5. The figure shows that the system
matrix is sparse, and banded; but as we have previously noted, matrices of this form do not admit a
sparse LU decomposition. Thus, in general, if direct elimination methods are employed, obtaining
solutions may require as much as O(N 3 ) arithmetic operations, where N = Nx Ny . In practical
problems N can easily be 106 , or even larger, and it is clear from this that iterative methods usually
must be employed. We shall first consider application of SOR, and we then present use of ADI in
a pseudo-time marching procedure.
5.4.1
Successive Overrelaxation
Recall from Chap. 1 that SOR is an accelerated form of GaussSeidel iteration, which in turn is
just a modification of Jacobi iteration. Because of the sparsity of our system matrix, it is better
to derive the SOR iteration equations specifically for (5.47), rather than to use the general SOR
formula of Chap. 1. Thus, for the (i, j)th equation of (5.47) solving for uij leads to the Jacobi
iteration equation
1 (m)
(m)
(m)
(m)
(m+1)
ui1,j + ui,j1 + ui,j+1 + ui+1,j h2 fij .
uij
=
4
It is interesting to observe that uij is computed as simply the average of its four nearest neighbors
whenever fij = 0.
Now we modify this iteration by always using the most recent iterate. From the ordering of the
(m+1)
(m+1)
(m+1)
uij shown in Fig. 5.5 we see that when uij
is to be computed, ui1,j and ui,j1 will already
be known. Thus, we arrive at the GaussSeidel iteration formula
1 (m+1)
(m+1)
(m)
(m)
(m+1)
ui1,j + ui,j1 + ui,j+1 + ui+1,j h2 fij .
uij
=
4
If we denote this result by uij , the SOR formula can be expressed as
(m+1)
uij
(m)
= uij + (1 )uij
(m)
(m)
(m)
uij + uij ,
= uij + uij uij
where is the relaxation parameter. (For problems of the type treated herein, 0 < < 2 must
hold.)
.
..
.
..
.
.
.
.
.
.
..
.
..
.
..
.
.
..
.
..
.
..
.
.
.
.
.
.
.
.
.
..
.
..
.
0
1
..
.
0
..
.
..
..
..
..
.
..
..
..
4
..
.
1
..
.
..
.
1
..
..
..
.
.
1
.
..
4
..
.
..
1
..
.
..
..
..
.
.
..
1
..
.
.
..
0
..
.
..
.
..
..
..
..
..
..
..
.
.
..
.
1
1
..
1
..
.
..
.
..
..
..
..
..
u11
g11 /h2
0
.
..
..
.
.
.
.
..
.. .
.
.
.
..
.
.
..
..
.
..
2
/h
g
u1,Ny
1,Ny
.
g /h2
..
u21
21
f
22
u22
..
. .
..
.
.
.. .
..
.
.. ..
ij
.. uij = h2
..
..
..
.
.
.. ..
.
..
..
..
.
.
g
/h
.. uNx ,1
Nx ,1
.
.
.
.
..
..
..
..
. .
.
.
..
.. .
.
.
..
.
.
.
0 .
2
1
gNx ,Ny /h
uNx ,Ny
146
147
It is well known that performance of SOR is quite sensitive to the value employed for the iteration
parameter . Optimal values of this parameter can be derived analytically for Poisson/Dirichlet
problems posed on rectangular domains (see, e.g., Young [39]). In fact, if the domain is the unit
square and the discretization employs uniform grid spacing h in both directions, it can be shown
that the optimal value of , denoted b , is given as
b =
1+
2
.
1 cos2 h
148
such difficulties are mitigated to some extent. Furthermore, the operation count per iteration for
SOR is far lower than that for essentially any other method, so although the convergence rate
of SOR may be inferior to that of other methods, it still is often used because its required total
arithmetic is often competitive, and it is easily coded.
5.4.2
Because of the convergence difficulties often experienced when SOR is applied to more realistic
problems, various alternatives have been studied. One of the early more successful of these was the
alternating-direction-implicit (ADI) algorithm. The treatment we shall present is very similar to
that given earlier for parabolic equations. In fact, we introduce a temporal operator, and formally
integrate the resulting time-dependent problem to steady state. Such an approach is sometimes
called the method of false transients. Its efficiency depends upon use of unconditionally stable
time-stepping procedures, so that relatively large time steps may be taken. While this would lead
to great inaccuracy in an actual transient problem, it is of no concern when only a steady solution
is required because all contributions to temporal truncation error vanish at steady state.
We begin again with Eqs. (5.45) and (5.46). But we now introduce a temporal operator corresponding to a pseudo time, , and write (5.45) as
u = uxx + uyy f (x, y) .
We already know that the PeacemanRachford method is unconditionally stable for linear problems
of this type. Thus, we can approximate the above with time-split formulas
k 2
k
k 2
(n)
(5.49a)
I D0,x ul,m = I + D0,y ul,m fl,m ,
2
2
2
k 2
k
k 2
(n+1)
(5.49b)
I D0,y ul,m = I + D0,x ul,m fl,m,
2
2
2
l, m, and use any desired value of time step, k. We observe that (5.49a,b) is identical to (5.40a,b)
except for the sign of f (which has been changed for notational consistency with (5.45)), and
our altered notation for the time step indices. The latter is due to viewing pseudo-time stepping
as being equivalent to iteration in the present case. In this regard, we also note that use of an
unconditionally stable time-stepping procedure implies convergence of the iteration procedure. In
particular, for the linear problems treated here, the connection between the amplification matrix of
the time-stepping algorithm and the Lipschitz constant of the iteration scheme is quite direct. As
a consequence, our only real concern with (5.49) is consistency. We wish to show that as n ,
(5.49) becomes a consistent approximation to (5.44).
To do this we note that from stability of the time stepping, and hence convergence of the
iteration procedure, we have
(n+1)
(n)
ul,m = ul,m = ul,m
as n , provided a steady solution to (5.45) exists. It then follows that (5.49) collapses to
2
2
D0,x
+ D0,y
ul,m = fl,m ,
which is the consistent difference approximation to (5.45). Details of proving this are left to the
interested reader.
Implementation of ADI for the steady case is only slightly different from that for transient
calculations. The first, and more obvious, difference is that a convergence test must now be incorporated into the algorithm, since convergence to steady state must be tested. Second, the effort
149
expended in coding temporary storage for time-accurate ADI is no longer needed; new iterates are
typically used as soon as they are available (a la GaussSeidel iteration) since temporal accuracy
is no longer an issue for steady BVP solutions.
Finally, the fact that the time-step parameter, r = k/2h2 , of the transient calculations is now
a somewhat arbitrary iteration parameter. It has been found that iterations converge much more
rapidly if this parameter is varied from one iteration to the next. In general, only a finite (and
usually fairly small) set of different iteration parameters is used; one application of each of these
on successive iterations until each has been used once is called a cycle. An algorithm implemented
so as to use a different iteration parameter on each iteration is said to constitute a nonstationary
iteration procedure. When only a finite set of parameters is used repeatedly, the method is said to
be cyclic. Thus, the ADI method just described is often termed cyclic ADI.
It should be noted that the cyclic ADI parameters for optimal performance are not easily
predicted for any but the simplest problems. Here we present a particular sequence applicable for
solving the Poisson equation on the unit square. We also point out that our definition of r is onehalf the usual definition. Thus, values of r obtained from the usual sequence, due to Wachspress
[7] must be divided by two before insertion into the preceding formulas. We first determine the
number of parameters making up the sequence. This is given as the smallest integer n0 such that
(n0 1)
21
.
tan
2(Nx 1)
Then the iteration parameter sequence is given by
n
n0 1
1
2
rn+1 =
cot
2 cos2 ( 2(Nx 1) )
2(Nx 1)
5.5
Hyperbolic Equations
We begin this section with a fairly complete, but elementary, treatment of the classical second-order
wave equation. We discuss some of the basic mathematics associated with this equation, provide
a second-order accurate difference approximation for solving it, and then analyze the consistency
and stability of this scheme. We then present an introduction to the first-order wave equation, and
extend this treatment to first-order hyperbolic systems. The reader is referred to Strikwerda [32]
for a comprehensive treatment of these topics.
5.5.1
(5.50)
150
where c is the wave speed. Here we will mainly consider the Cauchy problem for this equation.
That is, we let x (, ) and prescribe the initial conditions
u(x, 0) = f (x) ,
(5.51a)
ut (x, 0) = g(x) .
(5.51b)
For this case, (5.50) has the exact solution (known as the dLambert solution),
Z x+ct
1
g() d .
f (x + ct) + f (x ct) +
u(x, t) =
2
xct
(5.52)
An important property of this solution is that for any fixed and finite (x, t), the value of u is
determined from values of the initial data restricted to the finite interval [x ct, x + ct]. This
interval is called the domain of dependence of the point (x, t) and is depicted in Fig. 5.6.
The lines running between (x ct, 0) and (x, t), and (x + ct, 0) and (x, t), are respectively the
right-running and left-running characteristics through the point (x, t). These are shown extended
beyond (x, t) in the figure to indicate the region of influence of the point (x, t). The value of u
at any point in this shaded region depends on the value at (x, t), as well as elsewhere, of course,
and the solution at any point outside this region will be independent of u(x, t). The slope of the
characteristics is determined from the wave speed c as
1 1
.
= tan
c
( x,t )
Rightrunning
characteristic
Domain of
dependence
Region of
influence
Leftrunning
characteristic
xct
x+ct
Figure 5.6: Analytical domain of dependence for the point (x, t).
151
k2 n
n
n
n1
u
2u
+
u
um
.
m1
m
m+1
h2
2 n
2
n
n
n1
un+1
m = 2(1 )um + um+1 + um1 um .
(5.53)
We immediately observe that this is a three-level scheme. Hence, when n = 1, the left-hand side
of (5.53) will provide results at n = 2. Moreover, use of n = 1 leads to only one term containing
initial datathe last term on the right-hand side. Thus, we must have an independent means to
produce results for n = 1. This can be done using the initial data as follows. First, the n = 0 time
level is specified by the initial condition u(x, 0) = f (x). Now, from a Taylor expansion we have
u1m
u0m
+k
u
t
0
k2
+
2
m
2u
t2
0
+ O(k3 ) .
(5.54)
(5.55)
2u
t2
0
2u
x2
0
d2 f
dx2
1
(fm1 2fm + fm+1 ) + O(h2 ) .
h2
(5.56)
Substitution of (5.55) and (5.56) into (5.54) results in an equation for u1m in terms of prescribed
initial data f (x) and g(x):
1
u1m = (1 2 )fm + 2 (fm+1 + fm1 ) + kgm .
2
(5.57)
This provides the additional information needed to start calculations that are then continued using
(5.53).
We next discuss implementation of (5.53). It should be noted that in practice the initial data,
f (x) and g(x) will only be specified on a finite interval, rather than for x (, ). As Fig. 5.7
indicates, the length of this interval must be set by the amount of information (number of spatial
grid points in the solution) required at the final time. Clearly, as can be seen from the form of
the mesh star, one grid point is lost from each end of the grid at each time step just as occurred
for the Cauchy problem in the parabolic case discussed earlier. Thus, once k and h are set, the
number of grid points on the initial interval required to obtain a desired number on the final interval
can be determined. The algorithm in this case is quite similar to Algorithm 5.1, and we leave its
development as an exercise for the reader.
Consistency and Stability of Wave Equation Discretization
Now that we have considered calculation using the formula (5.53), it is worthwhile to analyze it
with respect to consistency and stability. To check consistency, we require the following Taylor
152
n+1
Mesh star for um
n+1
Numerical
characteristic
n
Numerical domain
of dependence
n1
m1
m+1
Figure 5.7: Numerical domain of dependence of the grid point (m, n + 1).
expansions:
k3
k4
k2
(utt )nm + (uttt )nm + (utttt )nm +
2
6
24
3
2
k
k4
k
= unm k(ut )nm + (utt )nm (uttt )nm + (utttt )nm
2
6
24
2
3
h
h
h4
= unm + h(ux )nm + (uxx )nm + (uxxx )nm + (uxxxx )nm +
2
6
24
3
4
2
h
h
h
= unm h(ux )nm + (uxx )nm (uxxx )nm + (uxxxx )nm
2
6
24
un+1
= unm + k(ut )nm +
m
n1
um
unm+1
unm1
k4
(utttt )nm + ,
12
h4
(uxxxx )nm + .
12
153
1
(utttt )nm k2 (uxxxx )nm h2 + = 0 .
12
(5.58)
Clearly, as h, k 0, the dominant truncation error goes to zero, and the differential equation is
recovered; hence, the scheme (5.53) is consistent with the wave equation. Moreover, the method is
second-order accurate in both space and time.
We analyze the stability of (5.53) via a von Neumann analysis. This might seem particularly
appropriate in the present case because there are no boundary conditions. On the other hand,
because (5.49) is a three-level difference scheme, it must be recognized that the von Neumann
condition supplies only a necessary (and not sufficient) stability requirement. The preferred way
to carry out analysis of a multi-level difference equation is to reduce it to a system of two-level
equations. For the present problem we define
n+1
vm
= unm .
(5.59a)
n+1
vm
= unm .
(5.59b)
and
n
+
+ unm
um+1 = I + h
x
2 x2
h2 2
F
7 I + h(i) + unm = eih unm ,
2
where F denotes Fourier transform. Then Eqs. (5.54) become
n
n+1
2
2 h
unm vm
um = 2 1 2 sin
2
n+1
vm
= unm ,
n+1
vm
1
h
2
1 unm
.
n
v
m
0
Now using arguments similar to those employed earlier when studying stability for difference approximations to parabolic equations, we know that the error vector must satisfy this same equation.
n R2 , we see that
Denoting this by zm
n+1
n
zm
= Czm
,
154
where
2 1 22 sin2
C
1
h
2
is the amplification matrix. As we noted earlier, the von Neumann necessary condition for stability
is
kCk 1 ,
where kk is the spectral norm. Hence, we need to calculate the eigenvalues of C. The characteristic
polynomial is
h
2
2 1 22 sin2
+1=0 ,
2
which has roots
1
2
h 2 2 h
h
.
2 sin
1
sin
= 1 2 sin
2
2
2
2
(5.60)
h
2
1 we have
(
1 )
2
h
h
h
i 2 sin
1 2 sin2
= 1 22 sin2
.
2
2
2
which implies |+ | = | | = 1. Second, if 2 sin2
instability. Hence, we must require
2 sin2
for stability. This implies that
2
h
2
h
1
2
1
sin2
h
2
and thus 2 1, 1. It follows that k h (or k h/c for c 6= 1) must hold for stability of
(5.54). We again emphasize that this is a necessary, but not sufficient, condition in this case.
The CFL Condition
We now discuss an aspect of hyperbolic equations that does not arise in either parabolic or elliptic
cases. Recall that we earlier, in Fig. 5.6, introduced the notion of domain of dependence of a point
value of the solution to a hyperbolic equation. Although we did not mention it previously, there
is an important relationship between the domain of dependence of the differential equation and
that of its difference approximation, shown in Fig. 5.7. These are not the same, in general. For
each domain of dependence, the respective solution is influenced only by the points included in
its domain. Thus, it follows that if the domain of dependence of the numerical approximation is
contained within the domain of dependence of the differential equation, initial data of the latter
could be arbitrarily changed outside the domain of dependence of the difference scheme (but inside
that of the differential equation), and the numerical solution would remain unchanged. In such a
155
case, the solution to the difference equation would not converge to that of the differential equation
even though it might be a consistent approximation. This is the content of the well-known CourantFriedrichs-Lewy (CFL) condition, often simply called the Courant condition.
Theorem 5.2 (CFL condition) In order for a difference approximation to a Cauchy problem for a
hyperbolic equation to converge (to the solution of the DE) as h, k 0, the domain of dependence
of the difference equation must include that of the differential equation.
For the wave equation we have been discussing, this implies that 1 must hold, which coincides
with the von Neumann condition.
Although we have not thus far treated implicit methods for hyperbolic problems, we note that
the CFL condition is satisfied, independent of the value of for such methods because the domain
of dependence of the difference equation is the entire spatial domain under consideration. We shall
not give any further treatment herein to the initial boundary value problem for the wave equation
except to note that implicit methods can be constructed quite readily merely by evaluating terms
from the spatial approximation at the advanced time. For example, we would have
1
1
n+1
n+1
n1
n
n+1
n+1
u
2u
+
u
u
2u
+
u
.
=
m
m
m
m
m1
m+1
k2
h2
It is easily seen that after rearrangement this leads to a tridiagonal system to be solved at each
time step. Initial conditions are treated just as in the explicit case discussed above, and boundary
conditions are handled in the usual way. We note, however, that implicit methods are not so widely
used in solving hyperbolic equations for two main reasons. First, the Courant condition for stability
is not overly restrictive, at least for 1-D problems (but in multidimensions stability requirements
for explicit schemes can become burdensome), and second, a CFL number ( = ck/h) close to unity
is generally required for accurate simulaton of time-dependent wave propagation in any case.
5.5.2
The last topic we shall treat here is one that is less classical than those discussed previously, but
one of great importance in practical computational physics: first-order hyperbolic equations. We
will first present several methods that are applicable to single equations, and then use one of these,
the LaxWendroff method, to solve a first-order hyperbolic system corresponding to the secondorder wave equation. Again, we will treat only the pure Cauchy problem; analysis of first-order
initial boundary value problems is well beyond the intended scope of these lectures. We refer the
interested reader to the monograph by Kreiss and Oliger [18] for a more thorough introduction to
this topic.
The First-Order Wave Equation
The first-order equation we shall study here is
ut + aux = 0 ,
a > 0 , constant ,
(5.61)
(5.62)
u(x, t) = f (x at) ,
(5.63)
156
as is easily checked. Associated with (5.63) is a single family of right-running characteristics with
slope
1 1
= tan
.
a
We require knowledge of this for construction of difference schemes which satisfy the CFL condition.
For the first-order equation (5.61), the CFL condition requires that the characteristics passing
through the advanced time level point pass inside the points of the mesh star at the current time
level. This is demonstrated in Fig. 5.8, and is equivalent to satisfaction of a CFL condition in the
sense of the CFL Theorem.
n +1
k
n
n 1
Characteristic
h
m 1
m +1
The simplest approximation to (5.61) satisfying the CFL condition is the first-order approximation
n
un unm1
un+1
m um
= a m
,
k
h
or
un+1
= unm a(unm unm1 ) = (1 a)unm + aunm1 ,
(5.64)
m
where = k/h. This scheme is explicit, and it is easily checked that the CFL condition is satisfied,
provided a 1. We defer the stability analysis to what will be developed below for another
method. An obvious disadvantage of this method is its first-order accuracy.
We now provide an approach that is second order in both space and time, known as the leap
frog scheme because it leaps over, i.e., does not include, the center point of its mesh star. It is
constructed by replacing derivatives in both space and time by centered differences. Thus, (5.61)
is approximated by
n1
un unm1
un+1
m um
+ a m+1
=0,
2k
2h
157
(5.65)
It is easily checked from the preceding figure that the CFL condition again requires a 1. The
stability analysis for (5.64) is identical in form to the von Neumann analysis earlier applied to
the second-order wave equation approximation; and the result is similar. Namely, von Neumann
stability coincides with satisfaction of the CFL condition. It is also easily shown by standard
methods (Taylor expansion) that the scheme is second-order accurate in both space and time.
The final method we shall consider for single first-order hyperbolic equations is the widely-used
in a Taylor series about unm to
LaxWendroff [20] scheme. To derive this method we expand un+1
m
second order, and then use the differential equation (5.65), to replace derivatives with respect to t
with those with respect to x. The derivation is completed by approximating the spatial derivatives
with centered differences. We have
un+1
= unm + k(ut )nm +
m
k2
(utt )nm + .
2
a2 k2
(uxx )nm + .
2
2 a2 n
a n
(um+1 unm1 ) +
(um1 2unm + unm+1 ) ,
2
2
or
a
a
(1 + a)unm1 + (1 2 a2 )unm (1 a)unm+1 .
(5.66)
2
2
Since this is only a two-level scheme, the stability analysis is analogous to that used earlier for
parabolic equations. We again note that the error in calculating the solution to (5.66) also satisfies
(5.65), and we represent this as
n
= n eimh .
zm
un+1
=
m
158
and
u2 = ux .
Then we have
(u1 )t (u2 )x = 0
from the original equation, and
(u2 )t (u1 )x = 0
from the transformation. These two equations can be represented as the first-order system
0 1 u1
1 0 u1
= 0,
(5.67)
+
1 0
u2 x
0 1 u2 t
or, equivalently,
Next set
u1
u2
0 1
A
1 0
0 1 u1
=0.
1 0 u2 x
t
and
U (u1 , u2 )T ,
to obtain
Ut AUx = 0 .
(5.68)
k2
(Utt )nm + .
2
1
n+1
n
Um
= Um
+ Ak(Ux )nm A2 k2 (Uxx )nm + .
2
Finally, we can write this in terms of centered-difference operators as follows:
k2 2 2
n
n+1
.
Um = I + kAD0 A D0 Um
2
(5.69)
This same procedure can be applied to various other hyperbolic systems with system matrices
that are different from the matrix A given here. Moreover, although it is not straightforward, the
LaxWendroff method can also be applied to variable-coefficient problems.
5.6. SUMMARY
5.6
159
Summary
In this chapter we have studied elementary numerical methods, based on the finite-difference approximation, for solving the partial differential equations of classical mathematical physics, as well
as for first-order wave equations and first-order hyperbolic systems.
For parabolic equations typified by the heat equation we presented both explicit and implicit
methods. We should recall that explicit methods are generally less stable than implicit ones, and
this is especially true in the parabolic case. Particularly for initial boundary value problems we
recommend use of the CrankNicolson method in most situations even though there is considerably
more required arithmetic per time step in comparison with explicit techniques because of the need to
solve tridiagonal linear systems at each time step. But far larger time steps can be used, so required
arithmetic is essentially always less than that for an explicit method. In the 2-D case we gave a
fairly detailed treatment of the PeacemanRachford ADI scheme applied to the heat equation.
This is a particular instance of a wide class of methods often referred to as time splitting. They
are extremely effective for solving time-dependent problems because of their favorable stability
properties, and because of the efficiency with which they can be solved, especially on modern
symmetric multi-processor (SMP) and cluster architectures.
We discussed only two very basic techniques for solving elliptic equations such as Laplaces
equation, namely SOR and ADI. We noted that there are by now many more sophisticated methods,
but the majority of these tend to be of somewhat limited applicability. Variants of SOR and ADI
are naturally highly parallelizable, and because of this they may yet prove to be as effective as
many of the modern procedures such as conjugate gradient, incomplete LU decompositions and
generalized minimum residual (GMRES) to name a few.
In any case, it is important to recognize that iterative methods (at least in the guise of pseudotime marching) are essentially always used for solving the sparse, banded linear systems arising
from discretization of elliptic PDEs. Such systems can be extremely large (106 106 , or larger,
matrices are now common), and direct methods requiring O(N 3 ) arithmetic are not usually viable
techniques. But it must be mentioned that direct methods exist that require only O(N 2 ) arithmetic
operations when applied to the sparse systems considered here; if these can be parallelized, then
they may prove to be effective alternatives to iteration in some cases.
Our discussions of hyperbolic PDEs included numerical treatment of the classical second-order
wave equation as well as first-order wave equations in scalar and vector form. Explicit methods
were emphasized because they are very efficient and easily implemented, and in the hyperbolic
case the requirements for stability (generally satisfaction of a Courant condition) tend to impose
far less restriction on the time step than is true for parabolic problems. We mainly treated only
the Cauchy problem. This was done to avoid the technicalities of numerical boundary conditions
for hyperbolic problems. This is still a very active area of research, particularly in the context of
computational electromagnetics, but also in computational fluid dynamics.
The main computational techniques we presented include straightforward centered differencing (in both space and time) for the second-order wave equation and the leap-frog method, the
equivalent scheme for first-order equations. In addition, we provided a fairly complete treatment
of the LaxWendroff scheme for first-order equations and systems. Although these are among the
most widely-used methods for hyperbolic equations, there are numerous other possibilities, as discussed in Strikwerda [32] and elsewhere. The interested reader is encouraged to consult the extant
literature.
160
References
[1] T. M. Apostol. Mathematical Analysis. Addison-Wesley Pub. Co., Reading, MA, second
edition, 1974.
[2] P. W. Berg and J. L. McGregor. Elementary Parial Differential Equations. Holden-Day, San
Francisco, 1966.
[3] G. Dahlquist. A Special Stability Problem for Linear Multistep Methods. BIT 3, New York,
1963.
[4] P. J. Davis and P. Rabinowitz. Methods of Numerical Integration. Academic Press, New York,
NY, 1975.
[5] C. deBoor. A Practical Guide to Splines. V. A. Barker (ed.), Springer-Verlag, New York, NY,
1978.
[6] C. B. Moler J. J. Dongarra, J. R. Bunch and G. W. Stewart. LINPACK Users Guide. SIAM,
Philadelphia, PA, 1979.
[7] J. Douglas, Jr. and J. E. Gunn. A general formulation of alternating direction methods, part
i. parabolic and hyperbolic problems. Numer. Math. 6, 428453, 1964.
[8] I. S. Duff, A. M. Erisman, and J. K. Reid. Direct Methods for Sparse Matrices. Clarendon
Press, Oxford, 1986.
[9] G. Forsythe and C. B. Moler. Computer Solution of Linear Algebraic Systems. Prentice-Hall,
Inc., Englewood Cliffs, NJ, 1967.
[10] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential Equations. PrenticeHall, Inc., Englewood Cliffs, NJ, 1971.
[11] K. E. Gustafson. Introduction to Partial Differential Equations and Hilbert Space Methods.
John Wiley and Sons, New York, 1980.
[12] W. Hackbusch. Iterative Solution of Large Sparse Systems of Equations. Springer-Verlag, NY,
1994.
[13] P. Henrici. Elements of Numerical Analysis. John Wiley & Sons, New York, 1964.
[14] R. W. Hornbeck. Numerical Methods. Quantum Publishers, Inc., New York, NY, 1975.
[15] E. Isaacson and H. B. Keller. Analysis of Numerical Methods. John Wiley & Sons, Inc., New
York, NY, 1966 (now available as the Dover paperback first published in 1994).
161
162
REFERENCES
[16] L. W. Johnson and R. D. Riess. Numerical Analysis. Addison-Wesley Pub. Co., Reading, MA,
second edition, 1982.
[17] H. B. Keller. Numerical Methods for Two-Point Boundary-Value Problems. Dover Pub., Inc.,
New York, NY, 1992.
[18] H. Kreiss and J. Oliger. Methods for the Approximate Solution of Time Dependent Problems.
GARP Pub, 1973.
[19] L. Lapidus and J. H. Seinfield. Numerical Solution of Ordinary Differential Equations. Academic Press, New York, NY, 1971.
[20] P. D. Lax and B. Wendroff. Systems of Conservation Laws, volume 13. Can. Pure Appl.
Math., 1960.
[21] A. R. Mitchell and D. F. Griffiths. The Finite Difference Method in Partial Differential Equations. John Wiley & Sons, Inc., Chichester, 1980.
[22] A. M. Ostrowski. Solution of Equations and Systems of Equations. Academic Press, New York,
second edition, 1966.
[23] D. W. Peaceman and H. H. Rachford, Jr. The numerical solution of Parabolic and Elliptic
Differential Equations, SIAM J. bf 3, 2841, 1955.
[24] R. D. Richtmyer and K. W. Morton. Difference Methods for Initial-Value Problems. John
Wiley & Sons, Inc., New York, NY, second edition, 1967.
[25] A. Ruhe. Computation of Eigenvalues and Eigenvectors in Sparse Matrix Techniques. V.
A. Barker (ed.), Springer-Verlag, Berlin, 1977.
[26] Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Co., Boston, 1996.
[27] B. T. Smith, P. Bjrstad and W. Gropp. Domain Decomposition, Parallel Multilevel Methods
for Elliptic Partial Differential Equations. Cambridge University Press, Cambridge, 1996.
[28] B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y. Ikebe, V. C. Klema and C.
B. Moler. Matrix Eigensystem Routines - EISPACK Guide, volume 6 of Lecture Notes in
Computer Science. Springer-Verlag, New York, NY, second edition, 1976.
[29] I. Stakgold. Boundary Value Problems of Mathematical Physics Volume I. Society for Industrial
and Applied Mathematics, Philadelphia, 2000.
[30] J. Stoer and R. Bulirsch. Introduction to Numerical Analysis. Springer-Verlag, New York,
1980.
[31] G. Strang and G. J. Fix. An Analysis of the Finite Element Method. Prentice-Hall, Inc.,
Englewood Cliffs, NJ, 1973.
[32] J. C. Strikwerda. Finite Difference Schemes and Partial Differential Equations. Wadsworth
& Brooks/Cole, Pacific Groove, CA, 1989.
[33] A. H. Stroud. Approximate Calculation of Multiple Integrals. Prentice-Hall, Inc., Englewood
Cliffs, NJ, 1971.
REFERENCES
163
[34] L. H. Thomas. Elliptic Problems in Linear Difference Equations over a Network, Watson
Sci. Comput. Lab. Rept., Columbia University, NY, 1949.
[35] J. F. Thompson, B. K. Soni and N. P. Weatherill (Eds.) Handbook of GRID GENERATION.
CRC Press, Boca Raton, 1999.
[36] J. F. Traub. Iterative Methods for the Solution of Equations. Prentice-Hall, Inc., Englewood
Cliffs, NJ, 1964.
[37] J. H. Wilkinson. The Algeibraic Eigenvalue Problem. University Press, London, 1965.
[38] J. H. Wilkinson and C. Reinsch. Linear Algebra. Springer-Verlag, New York, NY, 1971.
[39] D. M. Young. Iterative Solution of Large Linear Systems. Academic Press, New York, NY,
1971.