Numerical Analysis - Hardcover

Scott, Larkin Ridgway

Hardcover

ISBN 10: 0691146861 ISBN 13: 9780691146867

Publisher: Princeton University Press, 2011

Synopsis

Computational science is fundamentally changing how technological questions are addressed. The design of aircraft, automobiles, and even racing sailboats is now done by computational simulation. The mathematical foundation of this new approach is numerical analysis, which studies algorithms for computing expressions defined with real numbers. Emphasizing the theory behind the computation, this book provides a rigorous and self-contained introduction to numerical analysis and presents the advanced mathematics that underpin industrial software, including complete details that are missing from most textbooks.


Using an inquiry-based learning approach, Numerical Analysis is written in a narrative style, provides historical background, and includes many of the proofs and technical details in exercises. Students will be able to go beyond an elementary understanding of numerical simulation and develop deep insights into the foundations of the subject. They will no longer have to accept the mathematical gaps that exist in current textbooks. For example, both necessary and sufficient conditions for convergence of basic iterative methods are covered, and proofs are given in full generality, not just based on special cases.


The book is accessible to undergraduate mathematics majors as well as computational scientists wanting to learn the foundations of the subject.

• Presents the mathematical foundations of numerical analysis


• Explains the mathematical details behind simulation software


• Introduces many advanced concepts in modern analysis


• Self-contained and mathematically rigorous


• Contains problems and solutions in each chapter


• Excellent follow-up course to Principles of Mathematical Analysis by Rudin

About this title

About the Author

L. Ridgway Scott is the Louis Block Professor of Mathematics and Computer Science at the University of Chicago.

From the Back Cover

"Very few modern books can be compared with the present text as an introduction to the mathematical aspects of numerical analysis. This is a very interesting book that can be used not only as a textbook but also as a reference."--Doron Levy, University of Maryland

"This is a strong text, one that is both modern and provides historical perspective."--Benjamin Fearing Akers, University of Illinois at Chicago

Excerpt. © Reprinted by permission. All rights reserved.

Numerical Analysis

By L. Ridgway Scott

PRINCETON UNIVERSITY PRESS

Copyright © 2011 Princeton University Press
All right reserved.
ISBN: 978-0-691-14686-7

Contents

Preface.................................................................xiChapter 1. Numerical Algorithms.........................................1Chapter 2. Nonlinear Equations..........................................15Chapter 3. Linear Systems...............................................35Chapter 4. Direct Solvers...............................................51Chapter 5. Vector Spaces................................................65Chapter 6. Operators....................................................81Chapter 7. Nonlinear Systems............................................97Chapter 8. Iterative Methods............................................115Chapter 9. Conjugate Gradients..........................................133Chapter 10. Polynomial Interpolation....................................151Chapter 11. Chebyshev and Hermite Interpolation.........................167Chapter 12. Approximation Theory........................................183Chapter 13. Numerical Quadrature........................................203Chapter 14. Eigenvalue Problems.........................................225Chapter 15. Eigenvalue Algorithms.......................................241Chapter 16. Ordinary Differential Equations.............................257Chapter 17. Higher-order ODE Discretization Methods.....................275Chapter 18. Floating Point..............................................293Chapter 19. Notation....................................................309Bibliography............................................................311Index...................................................................323

Chapter One

Numerical Algorithms

The word "algorithm" derives from the name of the Persian mathematician (Abu Ja'far Muhammad ibn Musa) AlKhwarizmi who lived from about 790 CE to about 840 CE. He wrote a book, Hisab al-jabr w'al-muqabala, that also named the subject "algebra."

Numerical analysis is the subject which studies algorithms for computing expressions defined with real numbers. The square-root vy is an example of such an expression; we evaluate this today on a calculator or in a computer program as if it were as simple as y². It is numerical analysis that has made this possible, and we will study how this is done. But in doing so, we will see that the same approach applies broadly to include functions that cannot be named, and it even changes the nature of fundamental questions in mathematics, such as the impossibility of finding expressions for roots of order higher than 4.

There are two different phases to address in numerical analysis:

? the development of algorithms and

? the analysis of algorithms.

These are in principle independent activities, but in reality the development of an algorithm is often guided by the analysis of the algorithm, or of a simpler algorithm that computes the same thing or something similar.

There are three characteristics of algorithms using real numbers that are in conflict to some extent:

? the accuracy (or consistency) of the algorithm,

? the stability of the algorithm, and

? the effects of finite-precision arithmetic (a.k.a. round-off error).

The first of these just means that the algorithm approximates the desired quantity to any required accuracy under suitable restrictions. The second means that the behavior of the algorithm is continuous with respect to the parameters of the algorithm. The third topic is still not well understood at the most basic level, in the sense that there is not a well-established mathematical model for finite-precision arithmetic. Instead, we are forced to use crude upper bounds for the behavior of finite-precision arithmetic that often lead to overly pessimistic predictions about its effects in actual computations.

We will see that in trying to improve the accuracy or efficiency of a stable algorithm, one is often led to consider algorithms that turn out to be unstable and therefore of minimal (if any) value. These various aspects of numerical analysis are often intertwined, as ultimately we want an algorithm that we can analyze rigorously to ensure it is effective when using computer arithmetic.

The efficiency of an algorithm is a more complicated concept but is often the bottom line in choosing one algorithm over another. It can be related to all of the above characteristics, as well as to the complexity of the algorithm in terms of computational work or memory references required in its implementation.

Another central theme in numerical analysis is adaptivity. This means that the computational algorithm adapts itself to the data of the problem being solved as a way to improve efficiency and/or stability. Some adaptive algorithms are quite remarkable in their ability to elicit information automatically about a problem that is required for more efficient solution.

We begin with a problem from antiquity to illustrate each of these components of numerical analysis in an elementary context. We will not always disentangle the different issues, but we hope that the differing components will be evident.

1.1 FINDING ROOTS

People have been computing roots for millennia. Evidence exists [64] that the Babylonians, who used base-60 arithmetic, were able to approximate

v2 [approximately equal to] 1 + 24/60 + 51/60² + 10/60³ (1.1)

nearly 4000 years ago. By the time of Heron a method to compute square-roots was established [26] that we recognize now as the Newton-Raphson-Simpson method (see section 2.2.1) and takes the form of a repeated iteration

x <- 1/2 (x + y/x), (1.2)

where the backwards arrow <- means assignment in algorithms. That is, once the computation of the expression on the right-hand side of the arrow has been completed, a new value is assigned to the variable x. Once that assignment is completed, the computation on the right-hand side can be redone with the new x.

The algorithm (1.2) is an example of what is known as fixed-point iteration, in which one hopes to find a fixed point, that is, an x where the iteration quits changing. A fixed point is thus a point x where

x = 1/2 (x + y/x). (1.3)

More precisely, x is a fixed point x = f(x) of the function

f(x) = 1/2 (x + y/x), (1.4)

defined, say, for x [not equal to] 0. If we rearrange terms in (1.3), we find x = y/x, or x² = ^y. Thus a fixed point as defined in (1.3) is a solution of x² = y, so that x = ± vy].

To describe actual implementations of these algorithms, we choose the scripting syntax implemented in the system octave. As a programming language, this has some limitations, but its use is extremely widespread. In addition to the public domain implementation of octave, a commercial interpreter (which predates octave) called Matlab is available. However, all computations presented here were done in octave.

We can implement (1.2) in octave in two steps as follows. First, we define the function (1.4) via the code

function x=heron(x,y) x=.5*(x+y/x);

To use this function, you need to start with some initial guess, say, x = 1, which is written simply as

x=1

(Writing an expression with and without a semicolon at the end controls whether the interpreter prints the result or not.) But then you simply iterate:

x=heron(x,y)

until x (or the part you care about) quits changing. The results of doing so are given in table 1.1.

We can examine the accuracy by a simple code

function x=errheron(x,y) for i=1:5 x=heron(x,y); errheron=x-sqrt(y)

end

We show in table 1.1 the results of these computations in the case y = 2. This algorithm seems to "home in" on the solution. We will see that the accuracy doubles at each step.

1.1.1 Relative versus absolute error

We can require the accuracy of an algorithm to be based on the size of the answer. For example, we might want the approximation [??] of a root x to be small relative to the size of x:

[??]/x = 1 + d, (1.5)

where d satisfies some fixed tolerance, e.g., |d| = e. Such a requirement is in keeping with the model we will adopt for floating-point operations (see (1.31) and section 18.1).

We can examine the relative accuracy by the simple code

function x=relerrher(x,y) for i=1:6 x=heron(x,y); errheron=(x/sqrt(y))-1

end

We leave as exercise 1.2 comparison of the results produced by the above code relerrher with the absolute errors presented in table 1.1.

1.1.2 Scaling Heron's algorithm

Before we analyze how Heron's algorithm (1.2) works, let us enhance it by a prescaling. To begin with, we can suppose that the number y whose square root we seek lies in the interval [1/2;, 2]. If y < 1/2 or y > 2, then we make the transformation

[??] = 4^k y (1.6)

to get [??] [member of] [1/2 , 2], for some integer k. And of course v[??] = 2^k vy. By scaling y in this way, we limit the range of inputs that the algorithm must deal with.

In table 1.1, we showed the absolute error for approximating v2, and in exercise 1.2 the relative errors for approximating v2 and v1/2 are explored. It turns out that the maximum errors for the interval [1/2, 2] occur at the ends of the interval (exercise 1.3). Thus five iterations of Heron, preceded by the scaling (1.6), are sufficient to compute vy] to 16 decimal places.

Scaling provides a simple example of adaptivity for algorithms for finding roots. Without scaling, the global performance (section 1.2.2) would be quite different.

1.2 ANALYZING HERON'S ALGORITHM

As the name implies, a major objective of numerical analysis is to analyze the behavior of algorithms such as Heron's iteration (1.2). There are two questions one can ask in this regard. First, we may be interested in the local behavior of the algorithm assuming that we have a reasonable start near the desired root. We will see that this can be done quite completely, both in the case of Heron's iteration and in general for algorithms of this type (in chapter 2). Second, we may wonder about the global behavior of the algorithm, that is, how it will respond with arbitrary starting points. With the Heron algorithm we can give a fairly complete answer, but in general it is more complicated. Our point of view is that the global behavior is really a different subject, e.g., a study in dynamical systems. We will see that techniques like scaling (section 1.1.2) provide a basis to turn the local analysis into a convergence theory.

1.2.1 Local error analysis

Since Heron's iteration (1.2) is recursive in nature, it it natural to expect that the errors can be expressed recursively as well. We can write an algebraic expression for Heron's iteration (1.2) linking the error at one iteration to the error at the next. Thus define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.7)

and let e_n = x_n - x = x_n - vy. Then by (1.7) and (1.3),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.8)

If we are interested in the relative error,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.9)

then (1.8) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.10)

Thus we see that

the error at each step is proportional to the square of the error at the previous step;

for the relative error, the constant of proportionality tends rapidly to 1/2 . In (2.20), we will see that this same result can be derived by a general technique.

1.2.2 Global error analysis

In addition, (1.10) implies a limited type of global convergence property, at least for x_n > x = vy. In that case, (1.10) gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.11)

Thus the relative error is reduced by a factor smaller than 1/2 at each iteration, no matter how large the initial error may be. Unfortunately, this type of global convergence property does not hold for many algorithms. We can illustrate what can go wrong in the case of the Heron algorithm when x_n < x = vy.

Suppose for simplicity that y = 1, so that also x = 1, so that the relative error is [??]_n = e_n - 1, and therefore (1.10) implies that

[??]_n+1 = 1/2 (1 - x_n)²/x_n. (1.12)

As x_n -> 0, [??]_n+1 -> 8, even though |[??]_n| < 1. Therefore, convergence is not truly global for the Heron algorithm.

What happens if we start with x₀ near zero? We obtain x1 near . From then on, the iterations satisfy x_n > vy, so the iteration is ultimately convergent. But the number of iterations required to reduce the error below a fixed error tolerance can be arbitrarily large depending on how small x₀ is. By the same token, we cannot bound the number of required iterations for arbitrarily large x₀. Fortunately, we will see that it is possible to choose good starting values for Heron's method to avoid this potential bad behavior.

1.3 WHERE TO START

With any iterative algorithm, we have to start the iteration somewhere, and this choice can be an interesting problem in its own right. Just like the initial scaling described in section 1.1.2, this can affect the performance of the overall algorithm substantially.

For the Heron algorithm, there are various possibilities. The simplest is just to take x₀ = 1, in which case

[??]₀ = 1/x - 1 = 1/vy - 1. (1.13)

This gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.14)

We can use (1.14) as a formula for [??]₁ as a function of x (it is by definition a function of y = x²); then we see that

[??]₁(x) = [??]₁(1/x) (1.15)

by comparing the rightmost two terms in (1.14). Note that the maximum of [??]₁(x) on [2^-1/2, 2^1/2] occurs at the ends of the interval, and

[??]₁(v2) = 1/2 (v2 - 1)²/v2 = 3/4 v2 - 1 [approximately equal to] 0.060660 . (1.16)

Thus the simple starting value x₀ = 1 is remarkably effective. Nevertheless, let us see if we can do better.

1.3.1 Another start

Another idea to start the iteration is to make an approximation to the squareroot function given the fact that we always have y [member of] [1/2, 2] (section 1.1.2). Since this means that y is near 1, we can write y = 1 + t (i.e., t = y - 1), and we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.17)

Thus we get the approximation x [approximately equal to] 1/2 (y + 1) as a possible starting guess:

x₀ = 1/2 (y + 1). (1.18)

But this is the same as x₁ if we had started with x₀ = 1. Thus we have not really found anything new.

1.3.2 The best start

Our first attempt (1.18) based on a linear approximation to the square-root did not produce a new concept since it gives the same result as starting with a constant guess after one iteration. The approximation (1.18) corresponds to the tangent line of the graph of vy at y = 1, but this may not be the best affine approximation to a function on an interval. So let us ask the question, What is the best approximation to vy on the interval [1/2, 2] by a linear polynomial? This problem is a miniature of the questions we will address in chapter 12.

The general linear polynomial is of the form

f(y) = a + by. (1.19)

If we take x₀ = f(y), then the relative error [??]₀ = [??]₀(y) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.20)

Let us write e_ab(y) = [??]₀(y) to be precise. We seek a and b such that the maximum of |e_ab(y)| over y [member of] [1/2 , 2] is minimized.

Fortunately, the functions

e_ab(y) = a/vy + b vy - 1 (1.21)

have a simple structure. As always, it is helpful to compute the derivative:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.22)

Thus e'_ab(y) = 0 for y = a/b; further, e'_ab(y) > 0 for y > a/b, and e'_ab(y) < 0 for y < a/b. Therefore, e_ab has a minimum at y = a/b and is strictly increasing as we move away from that point in either direction. Thus we have proved that

min e_ab = min e_ab = e_ab(a/b) = 2 vab - 1. (1.23)

Thus the maximum values of |e_ab| on [1/2 , 2] will be at the ends of the interval or at y = a/b if a/b [member of] [1/2 , 2]. Moreover, the best value of e_ab(a/b) will be negative (exercise 1.10). Thus we consider the three values

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1.24)

(Continues...)

Excerpted from Numerical Analysisby L. Ridgway Scott Copyright © 2011 by Princeton University Press. Excerpted by permission of PRINCETON UNIVERSITY PRESS. All rights reserved. No part of this excerpt may be reproduced or reprinted without permission in writing from the publisher.
Excerpts are provided by Dial-A-Book Inc. solely for the personal use of visitors to this web site.

Edition details

• Publisher: Princeton University Press
• Publication date: 2011
• Language: English
• ISBN 10: 0691146861
• ISBN 13: 9780691146867
• Binding: Hardcover
• Number of pages: 344