COMPUTATIONAL
MATHEMATICS
Models, Methods, and Analysis
with MATLAB and MPI
© 2004 by Chapman & Hall/CRC
CHAPMAN & HALL/CRC
A CRC Press Company
Boca Raton London New York Washington, D.C.
COMPUTATIONAL
MATHEMATICS
ROBERT E. WHITE
Models, Methods, and Analysis
with MATLAB and MPI
© 2004 by Chapman & Hall/CRC
This book contains information obtained from authentic and highly regarded sources. Reprinted material
is quoted with permission, and sources are indicated. A wide variety of references are listed. Reasonable
efforts have been made to publish reliable data and information, but the author and the publisher cannot
assume responsibility for the validity of all materials or for the consequences of their use.
Neither this book nor any part may be reproduced or transmitted in any form or by any means, electronic
or mechanical, including photocopying, microfilming, and recording, or by any information storage or
retrieval system, without prior permission in writing from the publisher.
The consent of CRC Press LLC does not extend to copying for general distribution, for promotion, for
creating new works, or for resale. Specific permission must be obtained in writing from CRC Press LLC
for such copying.
Direct all inquiries to CRC Press LLC, 2000 N.W. Corporate Blvd., Boca Raton, Florida 33431.
Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are
used only for identification and explanation, without intent to infringe.
© 2004 by Chapman & Hall/CRC
No claim to original U.S. Government works
International Standard Book Number 1-58488-364-2
Library of Congress Card Number 2003055207

Preface xiii
Introduction xv
1 Discrete Time-Space Models 1
1.1 Newton Cooling Models . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Heat Di
usion in a Wire . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Di
usion in a Wire with Little Insulation . . . . . . . . . . . . . 17
1.4 Flow and Decay of a Pollutant in a Stream . . . . . . . . . . . . 25
1.5 Heat and Mass Transfer in Two Directions . . . . . . . . . . . . . 32
1.6 Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . 42
2 Steady State Discrete Mo dels 51
2.1 Steady State and Triangular Solves . . . . . . . . . . . . . . . . . 51
2.2 Heat Di
usion and Gauss Elimination . . . . . . . . . . . . . . . 59
2.3 Cooling Fin and Tridiagonal Matrices . . . . . . . . . . . . . . . 68
2.4 Schur Complement . . . . . . . . . . . . . . . . . . . . . . . . . . 77
2.5 Convergence to Steady State . . . . . . . . . . . . . . . . . . . . 86
2.6 Convergence to Continuous Model . . . . . . . . . . . . . . . . . 91
3 Poisson Equation Models 99
3.1 Steady State and Iterative Methods . . . . . . . . . . . . . . . . 99
3.2 Heat Transfer in 2D Fin and SOR . . . . . . . . . . . . . . . . . 107
3.3 Fluid Flow in a 2D Porous Medium . . . . . . . . . . . . . . . . . 116
3.4 Ideal Fluid Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3.5 Deformed Membrane and Steepest Descent . . . . . . . . . . . . 130
3.6 Conjugate Gradient Method . . . . . . . . . . . . . . . . . . . . . 138
v
© 2004 by Chapman & Hall/CRC
vi CONTENTS
4 Nonlinear and 3D Models 145

8 Classical Methods for Ax = d 313
8.1 Gauss Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . 313
8.2 Symmetric Positive Definite Matrices . . . . . . . . . . . . . . . . 318
8.3 Domain Decomposition and MPI . . . . . . . . . . . . . . . . . . 324
8.4 SOR and P-regular Splittings . . . . . . . . . . . . . . . . . . . . 328
8.5 SOR and MPI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333
8.6 Parallel ADI Schemes . . . . . . . . . . . . . . . . . . . . . . . . 339
9 Krylov Methods for Ax = d 345
9.1 Conjugate Gradient Method . . . . . . . . . . . . . . . . . . . . . 345
9.2 Preconditioners . . . . . . . . . . . . . . . . . . . . . . . . . . . . 350
9.3 PCG and MPI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 356
9.4 Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 360
9.5 GMRES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365
© 2004 by Chapman & Hall/CRC
CONTENTS vii
9.6 GMRES(m) and MPI . . . . . . . . . . . . . . . . . . . . . . . . 372
Bibliography 379
© 2004 by Chapman & Hall/CRC
List of Figures
1.1.1 Temperature versus Time . . . . . . . . . . . . . . . . . . . . . . 6
1.1.2 Steady State Temperature . . . . . . . . . . . . . . . . . . . . . . 7
1.1.3 Unstable Computation . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1 Di
usion in a Wire . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.2 Time-Space Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2.3 Temperature versus Time-Space . . . . . . . . . . . . . . . . . . . 15
1.2.4 Unstable Computation . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.5 Steady State Temperature . . . . . . . . . . . . . . . . . . . . . . 16
1.3.1 Di
usion in a Wire with csur = .0000 and .0005 . . . . . . . . . . 22

3.4.1 Ideal Flow About an Obstacle . . . . . . . . . . . . . . . . . . . . 123
3.4.2 Irrotational 2D Flow y
{
 x
|
= 0 . . . . . . . . . . . . . . . . . . 124
3.4.3 Flow Around an Obstacle . . . . . . . . . . . . . . . . . . . . . . 128
3.4.4 Two Paths to (x,y) . . . . . . . . . . . . . . . . . . . . . . . . . . 129
3.5.1 Steepest Descent norm(r) . . . . . . . . . . . . . . . . . . . . . . 137
3.6.1 Convergence for CG and PCG . . . . . . . . . . . . . . . . . . . . 144
4.2.1 Change in F
1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
4.2.2 Temperatures for Variable c . . . . . . . . . . . . . . . . . . . . . 158
4.4.1 Heat Diusion in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 167
4.4.2 Temperatures Inside a 3D Fin . . . . . . . . . . . . . . . . . . . . 170
4.5.1 Passive Solar Storage . . . . . . . . . . . . . . . . . . . . . . . . . 171
4.5.2 Slab is Gaining Heat . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.5.3 Slab is Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.6.1 Domain Decompostion in 3D . . . . . . . . . . . . . . . . . . . . 182
4.6.2 Domain Decomposition Matrix . . . . . . . . . . . . . . . . . . . 186
5.1.1 Infected and Susceptible versus Space . . . . . . . . . . . . . . . 196
5.2.1 Grid with Artificial Grid Points . . . . . . . . . . . . . . . . . . . 199
5.2.2 Infected and Susceptible at Time = 0.3 . . . . . . . . . . . . . . 203
5.3.1 Three Curves with Jumps . . . . . . . . . . . . . . . . . . . . . . 206
5.3.2 Restored 1D Image . . . . . . . . . . . . . . . . . . . . . . . . . . 213
5.4.1 Restored 2D Image . . . . . . . . . . . . . . . . . . . . . . . . . . 219
5.5.1 Value of American Put Option . . . . . . . . . . . . . . . . . . . 222
5.5.2 P(S,T-t) for Variable Times . . . . . . . . . . . . . . . . . . . . . 226
5.5.3 Option Values for Variable Volatilities . . . . . . . . . . . . . . . 226

6.3.3 Concentration at t = 17 . . . . . . . . . . . . . . . . . . . . . . . 256
6.4.1 Fan-out Communication . . . . . . . . . . . . . . . . . . . . . . 262
6.6.1 Space Grid with Four Subblocks . . . . . . . . . . . . . . . . . . 269
6.6.2 Send and Receive for Pro cessors . . . . . . . . . . . . . . . . . . 270
7.2.1 A Fan-in Communication . . . . . . . . . . . . . . . . . . . . . . 283
© 2004 by Chapman & Hall/CRC
List of Tables
1.6.1 Euler Errors at t = 10 . . . . . . . . . . . . . . . . . . . . . . . . 45
1.6.2 Errors for Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
1.6.3 Errors for Heat . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.6.1 Second Order Convergence . . . . . . . . . . . . . . . . . . . . . 96
3.1.1 Variable SOR Parameter . . . . . . . . . . . . . . . . . . . . . . . 104
3.2.1 Convergence and SOR Parameter . . . . . . . . . . . . . . . . . 113
4.1.1 Quadratic Convergence . . . . . . . . . . . . . . . . . . . . . . . . 149
4.1.2 Local Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 149
4.2.1 Newton’s Rapid Convergence . . . . . . . . . . . . . . . . . . . . 157
6.1.1 Truth Table for Bit Adder . . . . . . . . . . . . . . . . . . . . . . 239
6.1.2 Matrix-vector Computation Times . . . . . . . . . . . . . . . . . 243
6.2.1 Heat Di
usion Vector Times . . . . . . . . . . . . . . . . . . . . . 246
6.3.1 Speedup and E
!ciency . . . . . . . . . . . . . . . . . . . . . . . 252
6.3.2 HPF for 2D Di
usion . . . . . . . . . . . . . . . . . . . . . . . . 254
6.4.1 MPI Times for trapempi.f . . . . . . . . . . . . . . . . . . . . . . 262
6.5.1 Matrix-vector Product mflops . . . . . . . . . . . . . . . . . . . . 265
6.5.2 Matrix-matrix Product mflops . . . . . . . . . . . . . . . . . . . . 268
6.6.1 Processor Times for Di
usion . . . . . . . . . . . . . . . . . . . . 272
6.6.2 Processor Times for Pollutant . . . . . . . . . . . . . . . . . . . . 273

ATLAB is a registered trademark of The
MathWorks, Inc. For product information, please contact:
The MathWorks, Inc.
3 Apple Hill Drive
Natick, MA 01760-2098 USA
Tel: 508-647-7000
Fax: 508-647-7001
?A.
xiii
© 2004 by Chapman & Hall/CRC
Variations of Chapters 1-4 and 6 have been taught at North Carolina State
these activities c an be found at the Krell Institute, http://www.krellinst.org.
puting Center, http://www.ncsc.org, has provided the students with valuable
http://www.mathworks.com/Web: www.mathw orks.com
E-mail: [email protected]
xiv PREFACE
I thank my close friends who have listened to me talk about this e
ort, and
esp ecially Liz White who has endured the whole process with me.
Bob White, July 1, 2003
© 2004 by Chapman & Hall/CRC
Introduction
Computational science is a blend of applications, computations and mathemat-
ics. It is a mode of scientific investigation that supplements the traditional
lab oratory and theoretical methods of acquiring knowledge. This is done by
formulating mathematical models whose solutions are approximated by com-
puter simulations. By making a sequence of adjustments to the model and
subsequent computations one can gain some insights into the application area
under consideration. This text attempts to illustrate this process as a method
for scientific investigation. Each section of the first six chapters is motivated

chapters are revisited from a parallel implementation perspective.
xv
© 2004 by Chapman & Hall/CRC
codes; see [14] for
[17] for information about building your own multiprocessor via free "NPACI
xvi INTRODUCTION
lectures. Routine homework problems are assigned, and two projects are re-
work experiences. This forms a semester course on numerical modeling using
partial di
erential equations.
be studied after Chapter 1 so as to enable the student, early in the semester,
to become familiar with a high performance computing environment. Other
course possibilities include: a semester course with an emphasis on mathemat-
parallel computation using Chapters 1 and 6-9 or a year course using Chapters
1-9.
This text is not meant to replace traditional texts on numerical analysis,
matrix algebra and partial di
erential equations. It does develop topics in these
areas as is needed and also includes modeling and computation, and so there
is more breadth and less depth in these topics. One important component of
computational science is parameter identification and model validation, and this
requires a physical laboratory to take data from experiments. In this text model
assessments have been restricted to the variation of model parameters, mo del
evolution and mathematical analysis. More penetrating expertise in various
asp ects of computational science should be acquired in subsequent courses and
work experiences.
Related computational mathematics education material at the first and sec-
ond year undergraduate level can be found at the Shodor Education Founda-
tion, whose founder is Robert M. Pano
, web site [22] and in Zachary’s book

erent temperatures in
di
erent locations. Because heat flows from hot to cold regions, the subsequent
mo del will be more complicated. In Section 1.4 a similar model is considered,
and the application will be to the prediction of the pollutant concentration
in a stream resulting from a source of pollution up stream. Both of these
models are discrete versions of the continuous model that are partial dierential
equations. Section 1.5 indicates how these models can be extended to heat and
In the last section variations of the mean value theorem are used
to estimate the errors made by replacing the continuous model by a discrete
model. Additional introductory materials can be found in G. D. Smith [23],
and in R. L. Burden and J. D. Faires [4].
1.1 Newton Cooling Models
1.1.1 Introduction
Many quantities change as time progresses such as money in a savings account
or the temperature of a refreshing drink or any cooling mass. Here we will
be interested in making predictions about such changing quantities. A simple
mathematical model has the form
x
+
= dx + e where d and e are given real
numbers,
x is the present amount and x
+
is the next amount. This calculation is
usually repeated a number of times and is a simple example of an of algorithm.
A computer is used to do a large number calculations.
1
© 2004 by Chapman & Hall/CRC
mass transfer in two directions, which is discussed in more detail in Chapters

)10
h
=
This is a floating point numb er with base equal to 10 where {
1
is not equal
to zero,
{
l
are integers between 0 and 9, the exponent h is also a bounded
integer and g is an integer called the precision of the floating point system. As-
sociated with each real number,
{, and its floating point approximate number,
io({), is the floating point error, io({)  {. In general, this error decreases as
the precision,
g, increases. Each computer calculation has some floating point
or roundo
 error. Moreover, as additional calculations are done, there is an
accumulation of these roundo errors.
Example. Let
{ = 1=5378 and io({) = 0=154 10
1
where g = 3. The
roundo
 error is
io({)  { = =0022=
The error will accumulate with any further operations containing io({), for
example,
io({)
2

of cooling is
x
n+1
 x
n
= fk(x
vxu
 x
n
)
x
n+1
= (1  fk)x
n
+ fk x
vxu
= dx
n
+ e where d = 1  fk and e = fk x
vxu
=
© 2004 by Chapman & Hall/CRC
1.1. NEWTON COOLING MODELS 3
The long run solution should be the room temperature, that is,
x
n
should
converge to
x
vxu

modeled by
gx
n
+ e
x
n+1
 x
n
= gx
n
+ e
where g = d1= The model given in (1.1.1) is called a first order finite dierence
model for the sequence of numbers
x
n+1
. Later we will generalize this to a
sequence of column vectors where d will be replaced by a square matrix.
1.1.4 Method
The "iterative" calculation of (1.1.1) is the most common approach to solving
(1.1.1). For example, if
d =
1
2
> e = 2 and x
0
= 10, then
x
1
=
1

An alternative method is to use the following "telescoping" calculation and
the geometric summation. Recall the geometric summation
1 +
u + u
2
+ · · · + u
n
and (1 + u + u
2
+ · · · + u
n
)(1  u) = 1  u
n+1
© 2004 by Chapman & Hall/CRC
4 CHAPTER 1. DISCRETE TIME-SPACE MODELS
Or, for
u not equal to 1
(1 +
u + u
2
+ · · · + u
n
) = (1  u
n+1
)@(1  u)=
Consequently, if |u| ? 1, then
1 +
u + u
2
+ · · · + u

n2
+ d
2
e + de + e
.
.
.
=
d
n+1
x
0
+ e(d
n
+ · · · + d
2
+ d + 1)
= d
n+1
x
0
+ e(1  d
n+1
)@(1  d)
=
d
n+1
(x
0
 e@(1  d)) + e@(1  d)= (1.1.2)

x
4
 x = 4=375  4 = (
1
2
)
4
(10  4)=
1.1.5 Implementation
The reader should be familiar with the information in MATLAB’s tutorial. The
input segment of the M
ATLAB code fofdh.m is done in lines 1-12, the execution
is done in lines 16-19, and the output is done in line 20. In the following m-file
© 2004 by Chapman & Hall/CRC
1.1. NEWTON COOLING MODELS 5
w is the time array whose first entry is the initial time. The array | stores the
approximate temperature values whose first entry is the initial temperature.
The value of
f is based on a second observed temperature, |_revhu, at time
equal to k_revhu. The value of f is calculated in line 10. Once d and e have
been computed, the algorithm is executed by the for loop in lines 16-19. Since
the time step k = 1, q = 300 will give an approximation of the temperature
over the time interval from 0 to 300. If the time step were to be changed from 1
to 5, then we could change
q from 300 to 60 and still have an approximation of
the temperature over the same time interval. Within the for loop we could look
at the time and temperature arrays by omitting the semicolon at the end of the
lines 17 and 18. It is easier to examine the graph of approximate temperature
versus time, which is generated by the M
ATLAB command plot(t,y).

The first calculation is for this
f and k = 5 so that d = 1  fk = 60@65 and
e = fk @
the steady state room temperature, x
vxu
= 70.
The next calculation is for a larger
f = 2@13> which is computed from a new
second observed temperature of
x_revhu = 100 after k_revhu = 5 minutes.
In this case for larger time step
k = 10 so that d = 1  (2@13)10 = 7@13
and
e = fk70 = (2@13)10 70 = 1400@13. the
© 2004 by Chapman & Hall/CRC
70 = 350 65. Figure 1.1.1 indicates the expected monotonic decrease to
In Figure 1.1.2 notice that
6 CHAPTER 1. DISCRETE TIME-SPACE MODELS
Figure 1.1.1: Temperature versus Time
computed solution no longer is monotonic, but it does converge to the steady
state solution.
The model continues to degrade as the magnitude of
d increases. In the
This is consistent
with formula (1.1.2). Here we kept the same
f, but let the step size increase
to
k = 15 and in this case d = 1  (2@13)15 = 17@13 and e = fk70 =
(2
@13)1050 = 2100@13= The vertical axis has units multiplied by 10

n+1
= (1.1.3)
© 2004 by Chapman & Hall/CRC
Figure 1.1.3 the computed solution oscillates and blows up!
1.1. NEWTON COOLING MODELS 7
Figure 1.1.2: Steady State Temperature
Figure 1.1.3: Unstable Computation
© 2004 by Chapman & Hall/CRC
8 CHAPTER 1. DISCRETE TIME-SPACE MODELS
Next, we want to estimate the
dffxpxodwlrq huuru = X
n+1
 x
n+1
under the assumption that the roundo errors are uniformly bounded
|
U
n+1
|  U ? 4=
For ease of notation, we will assume the roundo errors associated with d and
e have been put into the U
n+1
so that X
n+1
= dX
n
+ e+ U
n+1
. Subtract (1.1.1)
and this variation of (1.1.3) to get

.
.
.
=
d
n+1
(X
0
 x
0
) + d
n
U
1
+ · · · + U
n+1
Now let u = |d| and U be the uniform bound on the roundo errors. Use the
geometric summation and the triangle inequality to get
|
X
n+1
 x
n+1
|  u
n+1
|X
0
 x
0
| + U(u

1.2. HEAT DIFFUSION IN A WIRE 9
minutes. Use the
f = 1/65 and k = 1.
6. We wish to calculate the amount of a savings plan for any month,
n> given a
fixed interest rate,
u, compounded monthly. Denote these quantities as follows:
x
n
is the amount in an account at month n, u equals the interest rate com-
pounded monthly, and g equals the monthly deposit. The amount at the end
of the next month will be the old amount plus the interest on the old amount
plus the deposit. In terms of the above variables this is with d = 1 + u@12 and
e = g
x
n+1
= x
n
+ x
n
u@12 + g
= dx
n
+ e=
(a). Use (1.1.2) to determine the amount in the account by depositing $100
each month in an account, which gets 12% compounded monthly, and over time
intervals of 30 and 40 years ( 360 and 480 months).
(b). Use a modified version of fofdh.m to calculate and graph the amounts
in the account from 0 to 40 years.
7. Show (1.1.5) follows from (1.1.4).

Compare the new curve with Figure
1.1.1.
10 CHAPTER 1. DISCRETE TIME-SPACE MODELS
change in time and
(change in temperature)/(change in space).
The last term is a good approximation provided the change in space is small,
and in this case one can use the derivative of the temperature with respect to
the single direction. The proportionality constant,
N, is called the thermal con-
ductivity. The N varies with the particular material and with the temperature.
Here we will assume the temperature varies over a smaller range so that
N is
approximately a constant. If there is more than one direction, then we must
replace the approximation of the derivative in one direction by the directional
derivative of the temperature normal to the surface.
Fourier Heat Law. Heat flows from hot to cold, and the amount of heat
transfer through a small surface area
D is proportional to the product of D> the
change in time and the directional derivative of the temperature in the direction
normal to the surface.
Consider a thin wire so that the most significant di
usion is in one direction,
{. The wire will have a current going through it s o that there is a source of
heat, i, which is from the electrical resistance of the wire. The i has units of
(heat)/(volume time). Assume the ends of the wire are kept at zero tempera-
ture, and the initial temperature is also zero. The goal is to be able to predict
the temperature inside the wire f or any future time and space location.
1.2.3 Model
In order to develop a model to do temperature prediction, we will discretize
both space and time and let

 x
n
l
)@k=
The heat diusing out the left face (when (x
n
l
 x
n
l
1
)@k A 0) is
D w N(x
n
l
 x
n
l
1
)@k.
Therefore, the heat from di
usion is
© 2004 by Chapman & Hall/CRC
This is depicted in the Figure 1.2.1 where the time step has not been indicated.
1.2. HEAT DIFFUSION IN A WIRE 11
Figure 1.2.1: Diusion in a Wire
D w N(x
n
l
+1

l
+1
 x
n
l
)@k  D w N(x
n
l
 x
n
l
1
)@k=
Now, divide by fDk, define  = (N@f)(w@k
2
) and explicitly solve for x
n+1
l
.
Explicit Finite Di
erence Model for Heat Diusion.
x
n+1
l
= (w@f)i + (x
n
l+1
+ x
n
l1

n+1
3
:
x
n+1
1
= (w@f)i + (x
n
2
+ 0) + (1  2)x
n
1
x
n+1
2
= (w@f)i + (x
n
3
+ x
n
1
) + (1  2)x
n
2
x
n+1
3
= (w@f)i + (0 + x
n
2