White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI

Подождите немного. Документ загружается.

Preface

This book evolved from the need to migrate computational science into under-

graduate education. It is intended for students who have had basic physics,

programming, matrices and multivariable calculus.

The choice of topics in the book has been inﬂuenced by the Undergraduate

Computational Engineering and Science Project (a United States Department

of Energy funded e

ort), which was a series of meetings during the 1990s.

These meetings focused on the nature and content for computational science

undergraduate education. They were attended by a diverse group of science

and engineering teachers and professionals, and the continuation of some of

University in fall semesters since 1992. The other four chapters were develop ed

in 2002 and taught in the 2002-03 academic year.

The department of mathematics at North Carolina State University has

given me the time to focus on the challenge of introducing computational science

materials into the undergraduate curriculum. The North Carolina Supercom-

tutorials and computer time on supercomputers. Many students have made

important suggestions, and Carol Cox Benzi contributed some course materials

with the initial use of M

AT LAB

. M

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

Variations of Chapters 1-4 and 6 have been taught at North Carolina State

these activities can 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.mathworks.com

E-mail: info@mathworks.com

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

Introduction

Computational science is a blend of applications, computations and mathemat-

ics. It is a mode of scientiﬁc investigation that supplements the traditional

lab oratory and theoretical methods of acquiring knowledge. This is done by

formulating mathematical mo dels 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 scientiﬁc investigation. Each section of the ﬁrst six chapters is motivated

by a particular application, discrete or continuous model, numerical method,

computer implementation and an assessment of what has been done.

Applications include heat di

usion to cooling ﬁns and solar energy storage,

pollutant transfer in streams and lakes, models of vector and multiprocessing

computers, ideal and porous ﬂuid ﬂows, deformed membranes, epidemic models

with dispersion, image restoration and value of American put option contracts.

The models are initially int roduced as discrete in time and space, and this allows

for an early introduction to partial di

erential equations. The discrete models

have the form of matrix products or linear and nonlinear systems. Methods in-

clude sparse matrix iteration with stability constraints, sparse matrix solutions

via variation on Gauss elimination, successive over-relaxation, conjugate gradi-

ent, and minimum residual methods. Picard and Newton methods are used to

approximate the solution to nonlinear systems.

Most sections in the ﬁrst ﬁve chapters have M

ATLAB

the very a

ordable current student version of MATLAB. They are intended

to be studied and not used as a "black box." The M

ATLAB codes should be

used as a ﬁrst step towards more sophisticated numerical modeling. These

codes do provide a learning by doing environment. The exercises at the end of

each section have three categories: routine computations, variation of models,

and mathematical analysis. The last four chapters focus on multiprocessing

algorithms, which are implemented using message passing interface, MPI; see

Rocks" cluster software. These chapters have elementary Fortran 9x codes

to illustrate the basic MPI subroutines, and the applications of the previous

chapters are revisited from a parallel implementation perspective.

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 imp ortant component of

computational science is parameter identiﬁcation 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, model

evolution and mathematical analysis. More penetrating expertise in various

asp ects of computational science should be acquired in subsequent courses and

work exp eriences.

Related computational mathematics education material at the ﬁrst 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

on programming [29]. Two general references for modeling are the undergradu-

ate mathematics journal [25] and Beltrami’s book on modeling for society and

biology [2]. Both of these have a variety of models, but often there are no

computer implemenations. So they are a good source of potential computing

projects. The book by Landau and Paez [13] has number of computational

physics models, which are at about the same level as this book. Slightly more

advanced numerical analysis references are by Fosdick, Jessup, Schauble and

Domik [7] and Heath [10].

The computer codes and updates for this book can be found at the web site:

The computer codes are mostly in M

ATLAB for Ch apters 1-5, and in Fortran

9x for most of the MPI codes in Chapters 6-9. The choice of Fortran 9x is

the author’s personal preference as the array operations are similar to those

in MATLAB. However, the above web site and the web site associated with

Pacheco’s book [21] do have C versions of these and related MPI codes. The

web site for this book is expected to evolve and also has links to sequences of

heat and pollution transfer images, book updates and new reference materials.

At North Carolina State University Chapters 1-4 are covered in 26 75-minute

quired, which can be chosen from topics in Chapters 1-5,related courses or

Chapter 6 on high performance computing can

ical analysis usingChapters 1-3, 8 and 9,asemester course with afocus on

http://www4.ncsu.edu/~white.

Chapter 1

Discrete Time-Space

Models

The ﬁrst three sections introduce diusion of heat in one direction. This is an

example of model evolution with the simplest model being for the temperature

of a well-stirred liquid where the temperature does not vary with space. The

model is then enhanced by allowing the mass to have di

erent temperatures in

erent locations. Because heat ﬂows 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 s ource 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

= 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.

mass transfer in two directions, which is discussed in more detail in Chapters

3 and 4.

2 CHAPTER 1. DISCRETE TIME-SPACE MODELS

Computers use a ﬁnite subset of the rational numbers (a ratio of two in-

tegers) to approximate any real number. This set of numbers may depend on

the computer being used. However, they do have the s ame general form and

are called ﬂoating point numbers. Any real number

{ can be represented by an

inﬁnite decimal expansion

{ = ±(={

· · · {

· · · )10

, and by truncating this we

can deﬁne the chopped ﬂoating point numbers.

Let

{ be any real number and denote a ﬂoating point number by

io({) = ±={

· · · {

= ±({

@10 + · · · + {

@10

)10

This is a ﬂoating point number with base equal to 10 where {

is not equal

to zero,

{

are integers between 0 and 9, the exponent h is also a bounded

integer and g is an integer called the precision of the ﬂoating point system. As-

sociated with each real number,

{, and its ﬂoating point approximate number,

io({), is the ﬂoating point error, io({)  {. In general, this error decreases as

the precision,

g, increases. Each computer calculation has some ﬂoating 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

where g = 3. The

roundo

 error is

io({)  { = =0022=

The error will accumulate with any further operations containing io({), for

example,

io({)

= =237 10

1

and

io({)

 {

= 2=37  2=36482884 = =00517116=

Repeated calculations using ﬂoating point numbers can accumulate signiﬁcant

roundo

 errors.

1.1.2 Applied Area

Consider the cooling of a well stirred liquid so that the temperature does not

dep end on space. Here we want to predict the temp erature of the liquid based on

some initial observations. Newton’s law of cooling is based on the observation

that for small changes of time,

k, the change in the temperature is nearly

equal to the product of the constant

f, the k and the dierence in the room

temperature and the present temperature of the co

ee. Consider the following

quantities: x

equals the temperature of a well stirred cup of coee at time w

vxu

equals the surrounding room temperature, and f measures the insulation

ability of the cup and is a positive constant. The discrete form of Newton’s law

of cooling is

n+1

 x

= fk(x

vxu

 x

)

n+1

= (1  fk)x

+ fk x

vxu

= dx

+ e where d = 1  fk and e = fk x

vxu

1.1. NEWTON COOLING MODELS 3

The long run solution should be the room temperature, that is,

should

converge to

vxu

as n increases. Moreover, when the room temperature is

constant, then x

should converge monotonically to the room temperature.

This does happen if we impose the constraint

? d = 1  fk>

called a stability condition, on the time step k. Since both f and k are positive,

d ? 1.

1.1.3 Model

The model in this case appears to be very simple. It consists of three constants

> d> e and the formula

n+1

= dx

+ e (1.1.1)

The formula must be used repeatedly, but with di

erent x

being put into the

right side. Often

d and e are derived from formulating how x

changes as n

increases (n reﬂects the time step). The change in the amount x

is often

modeled by

+ e

n+1

 x

= gx

+ e

where g = d1= The model given in (1.1.1) is called a ﬁrst order ﬁnite dierence

model for the sequence of numbers

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 =

> e = 2 and x

= 10, then

10 + 2 = 7

7 + 2 = 5

=5 + 2 = 4=75

=75 + 2 = 4=375

If one needs to compute

n+1

for large n, this can get a little tiresome. On

the other hand, if the calculations are being done with a computer, then the

ﬂoating point errors may generate signiﬁcant accumulation errors.

An alternative method is to use the following "telescoping" calculation and

the geometric summation. Recall the geometric summation

1 +

u + u

+ · · · + u

and (1 + u + u

+ · · · + u

)(1  u) = 1  u

n+1

4 CHAPTER 1. DISCRETE TIME-SPACE MODELS

Or, for

u not equal to 1

(1 +

u + u

+ · · · + u

) = (1  u

n+1

)@(1  u)=

Consequently, if |u| ? 1, then

1 +

u + u

+ · · · + u

+ · · · = 1@(1  u)

is a convergent geometric series.

In (1.1.1) we can compute

by decreasing n by 1 so that x

= dx

n1

+ e.

Put this into (1.1.1) and repeat the substitution to get

n+1

= d(dx

n1

+ e) + e

= d

n1

+ de + e

= d

(dx

n2

+ e) + de + e

= d

n2

+ d

e + de + e

n+1

+ e(d

+ · · · + d

+ d + 1)

= d

n+1

+ e(1  d

n+1

)@(1  d)

n+1

 e@(1  d)) + e@(1  d)= (1.1.2)

The error for the steady state solution, e@(1  d)> will be small if |d| is small,

or n is large, or the initial guess x

is close to the steady state s olution. A

generalization of this will be studied in Section 2.5.

Theorem 1.1.1 (Steady State Theorem) If d is not equal to 1, then the so-

lution of (1.1.1) has the form given in (1.1.2). Moreover, if |

d| ? 1, then the

solution of (1.1.1) will converge to the steady state solution

x = dx + e, that is,

x = e@(1  d). More precisely, the error is

n+1

 x = d

n+1

 e@(1  d))=

Example. Let d = 1@2> e = 2> x

= 10 and n = 3= Then (1.1.2) gives

3+1

= (1@2)

(10  2@(1  1@2)) + 2@(1  1@2) = 6@16 + 4 = 4=375=

The steady state solution is x = 2@(1 

) = 4 and the error for n = 3 is

 x = 4=375  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

ATLA B 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-ﬁle

1.1. NEWTON COOLING MODELS 5

w is the time array whose ﬁrst entry is the initial time. The array | stores the

approximate temperature values whose ﬁrst 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 t ime 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 temp erature

versus time, which is generated by the M

ATLAB command plot(t,y).

MATLAB Code fofdh.m

1. % This code is for the ﬁrst order ﬁnite dierence algorithm.

2. % It is applied to Newton’s law of cooling model.

3. clear;

4. t(1) = 0; % initial time

5. y(1) = 200.; % initial temperature

6. h = 1; % time step

7. n = 300; % number of time steps of length h

8. y_obser = 190; % observed temperature at time h_obser

9. h_obser = 5;

10. c = ((y_obser - y(1))/h_obser)/(70 - y(1))

11. a = 1 - c*h

12. b = c*h*70

13. %

14. % Execute the FOFD Algorithm

15. %

16. for k = 1:n

17. y(k+1) = a*y(k) + b;

18. t(k+1) = t(k) + h;

19. end

20. plot(t,y)

An application to heat transfer is as follows. Consider a cup of co

ee,

which is initially at 200 degrees and is in a room with temperature equal to

70, and after 5 minutes it cools to 190 degrees. By using k = k_revhu = 5,

= 200 and x

= x_revhu = 190, we compute from (1.1.1) that f = 1@65.

The ﬁrst 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 = 1 400@13. the

70 = 350 65. Figure 1.1.1 indicates the exp ected 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

k = 15 and in this case d = 1  (2@13)15 = 17@13 and e = fk70 =

@13)1050 = 2100@13= The vertical axis has units multiplied by 10

1.1.6 Assessment

Models of savings plans or loans are discrete in the sense that changes only occur

at the end of each month. In the case of the heat transfer problem, the formula

for the temperature at the next time step is only an approximation, which gets

better as the time step

k decreases. The co oling process is continuous because

the temperature changes at every instant in time. We have used a discrete

model of this, and it seems to give good predictions provided the time step is

suitably small. Moreover there are other modes of transferring heat such as

usion and radiation.

There may be signiﬁcant accumulation of roundo

 error. On a computer

(1.1.1) is done with ﬂoating point numbers, and at each step there is some new

roundo

 error

n+1

. Let X

= io(x

)> D = io(d) and E = io(e) so that

n+1

= DX

+ E + U

n+1

= (1.1.3)

Figure 1.1.3 the computed solution oscillates and blows up!