Elementary Theory and Application of Numerical Analysis: Revised Edition

By David G. Moursund, James E. Miller and Charles S. Duris

This updated introduction to modern numerical analysis is a complete revision of a classic text originally written in Fortran but now featuring the programming language C++. It focuses on a relatively small number of basic concepts and techniques. Many exercises appear throughout the text, most with solutions. An extensive tutorial explains how to solve problems with C++.

• Language: English
• Publisher: Dover Publications
• Release date: Apr 22, 2013
• ISBN: 9780486310398

Book preview

INDEX

PREFACE TO THE REVISED

2011 DOVER EDITION

This book is designed as a text for an undergraduate (junior-senior) level course in numerical analysis. A student that has an understanding of the calculus sequence has the necessary mathematics courses to succeed in this course. A course in differential equations and some computer programming will be helpful to the student but not essential. The book contains sufficient material for a one semester course. Some theory of the calculus is given, and the pure as well and the applied mathematics student can become further grounded in mathematical maturity as well as be challenged with algorithm development and computer programming. The language C++ is used to code the algorithms; however, the student can use the language of preference. Many of the problems in a course of this type lead to a great deal of computational effort. Because of this it will be to the student’s benefit to review his/her skills in computer programming early in the course. To help with this effort a short course in C++ is given in the appendix.

Much of the material in the text is basically algorithmic in nature. In general an algorithm is stated or derived and then illustrated by examples. Then the underlying theory is presented and discussed. Theorems giving conditions under which the algorithm is applicable are stated and in some cases proof of the theorems are given. Frequent numerical examples for applying the just presented algorithm are given. Where appropriate, C++ programs are included to show how the algorithm can be implemented on a computer. In general the treatment of topics from numerical analysis is consistent with that ordinarily used in a modern graduate-level course in numerical analysis. Thus, a firm foundation is laid for future studies in this area.

The foremost goal of this book is to provide an introduction to modern numerical analysis. A relatively small number of basic concepts and techniques are considered. The emphasis is upon how and why each of these methods works. Underlying the entire development is a consideration of the error-analysis aspects of the various problems and algorithms discussed.

A secondary goal of the text is to review and solidify some of the basic concepts from elementary calculus and to help raise the overall level of mathematical maturity of the student. Repeated use is made of the mean-value theorem, intermediate-value theorem and Taylor’s series. Considerable emphasis is placed upon theory and proofs, both in the text proper and in the exercises.

Throughout each chapter are found a large number of exercises. Answers to most of the exercises and short discussions on some of the theoretical questions are given at the end of the book.

A substantial portion of this text is derived from a 1967 book of the same name by David G. Moursund and Charles S. Duris. I used this text for many years in my numerical analysis course for junior and senior level students and found it to be sufficiently rigorous and able to cover the essential topics. When I found it no longer in print, I contacted David Moursund as Duris is now deceased. I was given permission to use the old text, however I thought best, and since this text is no longer in print and does such a good job of presenting the essential topics, I felt compelled to revise it and make it available. The programming in the text was with the language FORTRAN and has been replaced with the modern language C++.

I am grateful to my former colleague Dr. Tylene Garrett, Professor of Computer Science at Transylvania University for the technical review of the book. Finally, I am indebted to my wife Betty, for her understanding as I prepared the materials for this book.

James E. Miller

July, 2011

CHAPTER

1

Solution of Equations by

Fixed-Point Iteration

1-1 INTRODUCTION

It would seem appropriate to begin this text with a definition of numerical analysis. A useful, short definition of numerical analysis is as difficult to give as a useful, short definition of calculus. To know what calculus is, it is necessary to study the subject in some depth. The same is true of numerical analysis. In proceeding through this text, however, the reader might keep in mind the following definition due to Professor Preston C. Hammer: Numerical analysis is the effective representation of anything by anything. That is, the main goals of numerical analysis are to provide effective procedures for solving problems. Because many numerical procedures can be effectively implemented on a computer, computers play an important role in modern-day numerical analysis.

The main usefulness of a digital computer lies in its ability to follow a repetitive set of directions accurately and rapidly. Included in this capability is the ability to modify the basic set of directions being followed during the execution of the computer program. By an iterative procedure we mean a process whereby a basic set of directions is followed repeatedly; modification of these directions in the course of the execution is allowed. In this chapter, some basic iterative procedures for solving functional equations are developed. The discussion is limited to functions of one variable. The application of iterative procedures to the solution of a system of simultaneous equations is given in Chapter 2. Further application of iterative procedures to the numerical solution of differential equations is given in Chapter 7.

This chapter also seeks to review and solidify some of the basic results from elementary calculus. Continuity, limit, the mean-value theorem, and graphing are among the topics considered. The reader may find it helpful to refer to one of the standard texts in elementary calculus.

1-2 EXAMPLE

We shall use one example to illustrate most of the ideas and results of this chapter. Suppose that one wishes to solve the equation f(x) = 0, where

Here A > 0 is assumed to be a given constant. The peculiar form of f(x) will be explained in Sec. 1-6. Note that f(x) is not defined for x = 0. It can be seen that the equation has two solutions, x = ± A½. Thus the problem is essentially that of computing A½. The reader is undoubtedly familiar with one or more methods for accomplishing this. The procedure to be given here is easily programmed for a computer and is relatively fast.

Let g(x) = ½(x + A/x) so that the equation f(x) = 0 is equivalent to the equation x = g(x). Consider the following procedure. Let p0 > 0 be a first approximation to A½. Define p1 by p1 = g(p0). Define p2 by p2 = g(p1). In general,

Because the terms of the sequence {pn} are formed by successively substituting into g(x) the term most recently generated, Eq. (1-2) is known as a method of successive substitution. Let us apply this method to a particular example.

EXAMPLE 1-1. Suppose that A = 5 and po = 2. One measure of the accuracy of the approximation pn to 5½ is the quantity en = |pn² − 5|. A short table of {pn} and {en} is given below.

Notice the apparent rapid convergence of the sequence {en} to zero. In the course of this chapter we shall develop sufficient theory to prove that the sequence {pn} of Eq. (1-2) converges to A½.

EXERCISE 1. Suppose that A = 5 and p0 = 3 in Eq. (1-2). Compute p1, p2, p3, and e0, e1, e2, as in Example 1-1.

1-3 PRINCIPAL ALGORITHM

In this chapter we shall investigate one method for solving equations of the form x = g(x). Various modifications of this method, along with methods of writing equations of the form f(x) = 0 into the form x = g(x), will be considered.

Much of numerical analysis is devoted to the development and study of algorithms. By the term algorithm we shall mean a list of instructions specifying a sequence of operations to be used in solving a certain type of problem. Because the problems to be considered are mathematical, the operations used will be mathematical in nature.

An algorithm for approximating the square root of a positive number is given by Eq. (1-2). Ordinarily an algorithm will be accompanied by one or more theorems stating what problem the algorithm is designed to solve and when it works. The method of solving x = g(x) which we shall investigate is that of successive substitution, or fixed-point iteration. The terminology fixed point arises from the fact that one seeks a point p which remains fixed under application of the function g(x); that is, p = g(p). The method of fixed-point iteration is defined by Algorithm 1-1.

ALGORITHM 1-1. Let p0 be an approximation to a solution to the equation x = g(x). Generate the sequence {pn} recursively by the relation pn = g(pn−1), n = 1, 2….

FIG. 1-1 Flow chart for square-root algorithm.

There are a number of ways to state an algorithm. One way is in words and mathematical symbols, as above. Another way is by means of a flow chart. Still another way is by the use of a specially designed algorithmic language. The computer language ALGOL was designed explicitly for communication of algorithms. The reader is encouraged to code the algorithm in a language of choice and execute the code.

We shall give flow charts, and programs in some cases, and sample output for many of the algorithms discussed in this text. To illustrate, suppose that for arbitrary A, p0, and n one wishes to compute p1, p2, …, pn by the rule

In addition, the quantities |pi² − A| are to be computed. A flow chart is given in Fig. 1-1.

The flow chart shown in Fig. 1-1 provides for reading in the number A along with an initial approximation p and the number n of iterations to be performed. Printout is to include the data A, p, n and the elements of the sequences {pn} and {en} as in Example 1-1.

Program 1-1. C++ program for square-root algorithm

/* Sample Output

Input a, p and n:5 2 6

5 2 6

2 1

2.25 0.0625

2.23611 0.000193138

2.23607 1.46822e-007

2.23607 1.46822e-007

2.23607 1.46822e-007

Input a, p and n:10 3 5

10 3 5

3 1

3.16667 0.0277783

3.16228 1.98451e-005

3.16228 2.42537e-007

3.16228 2.42537e-007

*/

EXAMPLE 1-2. The quadratic formula may be thought of as an algorithm for solving equations of the form Ax² + Bx + C = 0. Suppose that it is desired to solve a sequence of quadratic equations where it is known that A ≠ 0 and B² − 4AC ≥ 0. The necessary square root is to be approximated by using 10 iterations of Eq. (1-2) with p0 = 1. (See Fig. 1-2.)

FIG. 1-2 Flow chart for solving a quadratic equation.

EXERCISE 2. Cramer’s rule may be thought of as an algorithm for solving simultaneous linear equations. State the algorithm for the case of a system of two simultaneous equations, and state a theorem which tells when the algorithm works.

EXERCISE 3. Let {x1, x2, …, xn} be a given set of real numbers (not necessarily distinct). State in flow-chart form an algorithm to perform each of the following problems.

(a) Determine the value of the largest element of the set.

(b) Determine the value of the smallest element of the set.

(c) Determine the subscript of the element of smallest subscript satisfying (a).

(d) Determine the subscript of the element of largest subscript satisfying (a).

EXERCISE 4. Let {x1, x2, …, xn} be a given set of real numbers. State in flow-chart form an algorithm which reads the set of numbers, orders them from largest to smallest, and prints out the numbers in that order.

EXERCISE 5. State in flow-chart form an algorithm for determining whether or not a given integer n > 1 is a prime number.

EXERCISE 6. State in flow-chart form an algorithm for determining both the mean and the median of a given set of N numbers.

EXERCISE 7. The algorithm given by the flow chart of Fig. 1-2 is not too useful as it stands. Add to it provisions to handle the case D = B²−4AC < 0. In addition, add a provision which stops the iteration when |P²−D| < 10−7 if this condition is satisfied before 10 iterations have been completed.

EXERCISE 8. Write a C++ program for the algorithm given in Exercise 7.

1-4 REVIEW OF BASIC CONCEPTS

Certain concepts from calculus arise repeatedly in numerical analysis. Two of these are limit and continuity. In this chapter we are interested in real-valued functions of a real variable and in real sequences. In Chapter 3 we shall be considering functions from Euclidean n space into Euclidean n space, where n > 1; in addition we shall discuss the possible convergence of sequences of points in Euclidean n space. Again in Chapter 7 we shall be faced with functions of more than one variable. Thus it behooves us to give precise definitions here and to state these definitions so that they easily extend to higher dimensions.

By definition, a function has a domain and a range. The domain is the set of points on which the function is defined. The range is the set of points which constitute the values of the function as it ranges over its entire domain. In discussing functions two other types of sets are frequently encountered.

DEFINITION 1-1. Let δ > 0 be arbitrary. The open interval (p − δ, p + δ) is called a neighborhood of the point p, with radius δ.

DEFINITION 1-2. Let δ > 0 be arbitrary. The set consisting of the two open intervals (p − δ, p) and (p, p + δ) is called a deleted neighborhood of the point p, with radius δ.

Algorithm 1-1 generates a sequence of points p1, p2, p3, …. In some cases this sequence will converge to a solution to the equation x = g(x), while in other cases it will diverge. Definition 1-3 of the limit of a sequence will be used in developing the theory of fixed-point iterative procedures.

DEFINITION 1-3. We write

to mean that, for each ε > 0, there is an integer N (which often depends on ε) such that | pn − p| < ε whenever n > N.

EXAMPLE 1.3. Suppose that pn = n/(2n + 3). Then

To prove this, suppose that ε > 0 is arbitrary. We consider

for large n. This can be written as

Observe that

for positive n. It is now a simple matter to select N such that n > N implies that

Specifically, if N ≥ 3/(4ε) is an integer and n > N, then

EXERCISE 9. Prove from the definition that

EXERCISE 10. Prove from the definition that

Next we define the limit of a function at a point.

DEFINITION 1-4. Suppose that f(x) is defined on a domain D. Let p be a point such that every deleted neighborhood of p contains points of D. Then we write

to mean that for each ε > 0 there exists a δ > 0 (which often depends on ε) such that | f(x) − v| < ε whenever 0 < | x − p | < δ and x D.

Notice in Definition 1-4 that we are concerned only with points in a deleted neighborhood of the point p. The point p need not be in the domain of definition of the function f. If p ∈ D, then we can inquire into what is called the continuity of f(x) at p.

DEFINITION 1-5. Suppose that f(x) is defined on a domain D and p ∈ D. We say that f(x) is continuous at p if

Definition 1-5 actually requires three things. In order to be continuous at a point p, the function f(x) must be defined at p. Next, the function f(x) must be such that

exists. Finally, the value of the above limit must be f(p).

EXAMPLE 1-4. Using the appropriate definitions, we shall prove that f(x) = 2x + 5 is continuous at x = 3.

First we note that f(3) =11. We must prove that

exists and has the value 11. Let ε > 0 be arbitrary. We examine the quantity

|(2x + 5) − 11|

for x near 3. This quantity may be simplified and written in the form

|(2x + 5) − 11| = |2x − 6| = 2|x − 3|

This is to be made less than ε by making |x − 3| less than some δ. It is now evident that if we set δ = ε/2 then |x − 3| < δ implies that 2|x − 3| < ε.

EXERCISE 11.

(a) Let f(x) = 3x [that is, prove that f(x) is continuous at x = −1],

(b) Prove that the above function is continuous for all x.

EXERCISE 12.

(a) Prove that the function f(x) = x² + 1 is continuous at the point x = 2.

(b) Prove that the above function is continuous for all x ≥ 1.

The reader may be familiar with alternative definitions of continuity. Typical is the ε, δ definition.

EXAMPLE 1-5 Suppose that f(x) is defined on D and is continuous at a point p ∈ D. Let ε > 0 be arbitrary. Show that there exists a δ > 0 such that

|f(x) − f(p)| < ε

whenever |x − p| < δ and x D.

Solution From Definition 1-5 we know that

Using the definition of limit, for the above ε > 0, we know there exists a δ > 0 such that

|f(x) − f(p)| < ε

whenever |x − p| < δ and x ∈ D. This completes the proof.

Frequently one uses a hypothesis that a function is continuous at each point of a given set. If f(x) is continuous at each point of a set S, we say that f(x) is continuous on S and write

f(x) ∈ C(S)

That is, C(S) designates the class of continuous functions on S, and f(x) belongs to C(S). If the nth derivative of f(x) is continuous on a set S, for some n ≥ 1, we write

f(x) ∈ Cn(S)

Cn(S) is the class of functions with n continuous derivatives on S.

Next, let us recall an important theorem from elementary calculus. The majority of functions we shall encounter are differentiable.

DEFINITION 1-6. Suppose that f(x) is defined on a domain D and p ∈ D.

Then f(x) is differentiable at p if

exists. The value of the limit, if it exists, is denoted by f′(p).

THEOREM 1-4-1. If f(x) is differentiable at a point p, it is continuous at p. The proof of Theorem 1-4-1 is left as an exercise. The reader may find it helpful to refer to an elementary calculus text to refresh his memory on the details of a proof of this important result.

EXERCISE 13. Prove Theorem 1-4-1.

EXERCISE 14.

(a) Prove that f(x) = x½ is continuous on [1,3].

(b) Is f(x) = x½ continuous on [0,1]? Justify your assertion. (Hint: Consider the point x = 0 separately.)

EXERCISE 15. Let R denote the real line, and let

f(x) = a0 + a1x + a2x² +…+ amxm

Prove that f(x) ∈ C(R) and f(x) ∈ Cn(R) for n = 1, 2, ….

EXERCISE 16. Suppose that f(x) ∈ C¹(S) and g(x) ∈ C¹(S). Prove that f(x) • g(x) ∈ C¹(S).

We now have enough definitions to prove an important result. Algorithm 1-1 can be applied to any equation written in the form x = g(x). Under certain conditions the sequence {pi} generated will converge. In such a case Theorem 1-4-2 is applicable.

THEOREM 1-4-2. and g(x) is continuous at the point p. Then p = g(p).

Proof. Theorem 1-4-2 says that, if the sequence generated by the algorithm converges and the function g(x) is continuous at the limit point, then the limit point is a solution to the equation x = g(x).

Let ε > 0 be arbitrary. From the definition of continuity, applied at the point p, there exists a δ > 0 such that if | x − p | < δ then |g(x) − g(p)| < ε. Next, apply the definition of limit to the sequence {pn}. For the δ given above there exists an N such that n > N implies that |pn − p| < δ. Thus, if n > N, it follows that | g(pn) − g(p)| < ε. We have now proved that

To complete the proof, note that g(pn)=pn+1. Therefore,

, the proof is complete.

The intermediate-value theorem states that if a function f is continuous on a closed interval [a, b], then f takes on all values between f(a) and f(b) as x ranges over [a, b]. An equivalent form of this theorem is given below. This theorem is often useful in proving that an equation has a solution.

THEOREM 1-4-3. Let f(xC([a, b]), and suppose that f(a) • f(b) < 0. That is, f(a) and f(b) are of opposite sign. Then there exists a point p (a, b) such that f(p) = 0.

EXAMPLE 1-6. Prove that the equation x − cos(x) = 0 has a solution.

Solution. Let f(x) = x − cos(x). Because f′(x) = 1 + sin(x) exists for all x, f(x) is continuous for all x. We shall use Theorem 1-4-3. It is easy to see that f(0) = −1 and f(π/2) = π/2. Hence f(x) = 0 has a solution in (0, π/2).

We could arrive at the same conclusion by considering the graph of f(x). However, it is easier to consider two separate graphs

y = x     y = cos x

Each point common to these two curves is a solution to the equation x = cos(x). Conversely, each solution to x = cos(x) is a point of intersection of the two curves. The graphs are given in Fig. 1-3.

FIG. 1-3

EXERCISE 17. Determine which of the equations below have at least one solution. For those which do, find an interval containing one of the solutions.

(a) x − e −x = 0

(b) x − ex = 0

(c) sin x − 2 cos x = 0

(d) x − a−x = 0, for fixed a ≥ 1

EXERCISE 18. Let P(x) = a0x³ + a1x² + a2x + a3, where the