N u m e r i c al M e t h o d s, S o f t w a r e, a nd A n a l y s is Second Edition John R. Rice Department of Computer Science Purdue University West Lafayette, Indiana ACADEMIC PRESS, INC. Harcourt Brace Jovanovich, Publishers Boston San Diego New York London Sydney Tokyo Toronto This book is printed on acid-free paper. © Copyright © 1993, 1983 by Academic Press, Inc. All rights reserved. No part of this publication may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopy, recording, or any information storage and retrieval system, without permission in writing from the publisher. IMSL and PROTRAN are the trademarks of IMSL Inc. for its proprietary computer software. No materials describing such software may be produced or distributed without the written permission of the owners of the trademark and license rights in the software and the copyright in the published materials. ACADEMIC PRESS, INC. 1250 Sixth Avenue, San Diego, CA 92101-4311 United Kingdom edition published by ACADEMIC PRESS LIMITED 24-28 Oval Road, London NW1 7DX ISBN 0-12-587755-2 Printed in the United States of America 92 93 94 95 EB 9 8 7 6 5 4 3 2 1 PREFACE Objectives Modem science depends on our capability to carry out substantial computations to explore the consequence of the laws of nature. These laws are expressed as mathematical models, and numerical methods are used throughout the computation. The objective of this text is to introduce science and engineering students to the methods, tools, and ideas of numerical computation. These students are assumed to have taken the standard introductory calculus sequence, plus an elementary introduction to matrix theory, ordinary differential equations, and Fortran programming. An introductory course in numerical methods faces a fundamental problem: students who want to do realistic computations need to know more than they can learn in one course. This text has about twice as much material as can be covered in a typical one-term course, and yet most topics are not covered in depth. The use of high-quality mathematical software is the leading way out of the dilemma of too little time to learn too much knowledge and know-how. The mathematical software approach is based on the premise that experts in numerical computation can incorporate enough of their experience and know-how into a computer program so that it performs better than the program a typical scientist can create with a reasonable effort. The validity of this approach has been established in many areas of numerical computation; the objective of this text may be restated as to present scientific problem solving using standard mathematical software. Organization and Content This book is organized into three parts, as shown in the diagram below. Part I is background and preliminary material. Few instructors will cover this thoroughly; in a one-term course it is feasible to summarize the main ideas in a few lectures. Part II is the core of the course, presenting the principal methods and ideas of numerical computation. Part III contains very relevant, but somewhat independent, material about software engineering and performance evaluation. There is rarely time in a one-term course to cover this material in detail, yet it is important that young scientists gain some appreciation for this material. xi xii PREFACE PART I: Preliminaries 1. MATH AND CS BACKGROUND 2. NUMERICAL SOFTWARE 3. ERRORS PART II: Core Material PART III: Software Engineering 12. SOFTWARE PRACTICE 4. MODELS AND FORMULAS 13. PERFORMANCE EVALUATION 5. INTERPOLATION 14. SOFTWARE VALIDATION 6. MATRIX COMPUTATION 15. PROTRAN 7. DERIVATIVES/INTEGRALS 8. NONLINEAR EQUATIONS 9. DIFFERENTIAL EQUATIONS 10. PARTIAL DIFFERENTIAL EQUATIONS 11. APPROXIMATION A fairly uniform approach is used in each area of numerical computation. First, an intuitive development is made of the problems and the basic methods for their solution. Then, relevant mathematical software is reviewed and its use outlined. In many areas rather extensive examples and case studies are given. Finally, a deeper analysis of the methods is presented as in traditional numerical analysis texts. The use of mathematical software tools may resolve the dilemma facing the practicing scientist, but it still leaves a dilemma for the instructor using this text: should one emphasize how to use the software or understanding how it works? This is not an easy choice to make; ample material is presented so that an instructor may choose whatever balance seems appropriate for the students' needs and the instructor's taste. Some sections and chapters are "starred" to help organize the material. The double-starred items contain material that is either more advanced (for example, osculatory interpolation and adaptive quadrature) or somewhat peripheral to the main topics (for example, sources of matrix problems and reporting computational experiments). Since the remaining material is still more than can be covered in a one-term course, the single-starred items are suggested as candidates to be skipped. Some single-starred chapters also have single-starred sections; this provides suggested topics for an instructor who wants to cover these chapters briefly. The dependence of the chapters is shown in the following diagram. PREFACE xiii Software This text discusses numerous programs and software packages, but it is not intended to be a user's guide for this software. The IMSL library and ACM Algorithms are discussed systematically as they are the most widely available general collections of high-quality mathematical software. It is not essential to have this particular software; any reasonably complete program library will have programs similar to those used in most examples and problems. We note (see Chapter 12) that the cost of developing just one good program is more than the cost of acquiring several hundred good programs from IMSL, ACM, and similar sources. These sources should be viewed as providing the initial step toward a software parts technology which might dramatically reduce the present heavy burden of software development costs. This text also uses the PROTRAN system introduced by IMSL. PROTRAN provides mathematical problem-solving statements which substantially ease the use of library routines and simulta- neously improve programming reliability. PROTRAN has two important additional attractions for use in this text. First, example programs in PROTRAN are much shorter than in Fortran, and even if PROTRAN is not available to the students, they can easily visualize how an equivalent Fortran program would be. Second, PROTRAN is an example of the very high level languages that most students will be using in a few years as this methodology spreads throughout scientific computation. Acknowledgments The most important acknowledgment is to the students at Purdue University who use preliminary versions of this text. Their suggestions and patience with typos and rough spots are greatly appre- ciated. Second, I thank Thomas Aird especially, for carefully reviewing the discussion of IMSL software and for making the PROTRAN system available. I thank Robert Lynch for his invaluable xiv PREFACE assistance in developing the library applications presented here; I also thank Dinesh Manocha and Thomas Stiller for their help in implementing these applications. I thank Carl de Boor for permission to include the text of two programs (INTERV and Example 5.6) from PPPACK and for many valuable insights into the nature of scientific computation. I thank John Reid for material used in the cover design and Figure 6.3. Many valuable comments were received from Neil Berman, Ronald Boisvert, Henry Darilek, Wayne Dyksen, Don Kainer, Iftekhar Karimi, Robert Lynch, and Granville Sewell. I thank them and the others who have contributed to the development of this text. I thank Doris Moore for typing the text and several revisions, Richard Gaudette for preparing the art work, and Gary Lelvis for coordinating the preparation, typing, and typesetting of the text. John R. Rice June, 1992 **1 MATHEMATICS AND COMPUTER SCIENCE BACKGROUND The background assumed in this text is a basic sequence in calculus, an introduction to differential equations and linear algebra (or matrices), and a working knowledge of Fortran programming. The purpose of this chapter is to give a highly condensed review of this background. The material provides references for the rest of the text, and the problems allow instructors or students to check their knowledge against the assumed background in calculus and linear algebra. Many instructors will choose to spend little, if any, time presenting this material in lectures. Some material in this book requires more background and experience than presented here, for example, the chapter on partial differential equations. An instructor must judge whether the students have the necessary experience and maturity for such material; some instructors might choose to use the material here as a first introduction to the topic. The sections marked by ** contain material which is more advanced or peripheral and which can be skipped; all the material requiring more background is in such sections. 1.1 CALCULUS The principal topics in calculus are the real and complex number systems, the concept of limits and convergence, and the properties of functions. Convergence of a sequence of numbers x, is defined as follows: The sequence X{ converges to the limit x* if, given any tolerance c > 0, there is an index Ν = N(c) so that for all i> Ν we have |a?- - χ* | < e. The notation for this is t lim X{ = x* Ϊ—+ΌΟ Convergence is also a principal topic of numerical computation, but with a different emphasis. In calculus one studies limits and convergence with analytic tools; one tries to obtain the limit or to show that convergence takes place. In computations, one has the same problem but little or no theoretical 1 2 MATHEMATICS AND COMPUTER SCIENCE BACKGROUND knowledge about the sequence. One is frequently reduced to using empirical or intuitive tests for convergence; often the principal task is to actually estimate the value of the tolerance c for a given The study of functions in calculus revolves about continuity, derivatives, and integrals. A function f(x) is continuous if lim./(*,•) = /(**) holds for all x* and all ways for the Xi to converge to x*. We list six theorems from calculus which are useful for estimating values that appear in numerical computation. Theorem 1.1 (Mean value theorem for continuous functions) Let f(x) be continuous on the interval [a, b]. Consider points XHI and XLOW in [a, b] and a value y so that f(XLOW) < y < f(XHI). Then there is a point in [a, 6] so that Theorem 1.2 (Mean value theorem for sums) Let f(x) be continuous on the interval [a, b], let x\,X2^..,xbe points in [a, 6], and let w\,W2,.. - ,w be positive numbers. Then there is a point n n £ in [a, b] so that η η £>*/(**) =/(Ο 5>i Theorem 1.3 (Mean value theorem for integrals) Let f(x) be continuous on the interval [a, b] and let w(x) be a nonnegative function [w(x) > 0] on [a, 6]. Then there is a point ξ in [a, 6] so that / w(x)f(x)dx = /(£) / w(x)dx Ja J a Theorems 1.2 and 1.3 show the analogy that exists between sums and integrals. This fact derives from the definition of the integral as / f(x)dx= lim y2f(xi)(xM ~Xi) Ja max\xi-xi\-+0 i+ where the points a?», with X( < x \, are a partition of [a, 6]. This analogy shows up for many nu- i+ merical methods where one variation applies to sums and another applies to integrals. Theorem 1.2 is proved from Theorem 1.1, and then Theorem 1.3 is proved by a similar method. The assumption that w(x) > 0(wi > 0) may be replaced by w(x) < 0(w{ < 0) in these theorems; it is essential that w(x) be of one sign as shown by the example w(x) = f(x) = χ and [a, b] = [-1,1]. Theorem 1.4 (Continuous functions assume max/min values) Let f(x) be continuous on [a, b] with |a|, |6| < oo. Then there are points XHI and XLOW in [a, b] so that for all χ in [a, 6] f(XHI) < f(x) < f(XLOW) The derivative of f(x) is defined by £ , / ( «+ *>-/(«> =/ ( x ) = U m dx Λ -Ο h CALCULUS 3 As an illustration of the difference between theory and practice, the quantity [/(x + ft) - f(x)]/h can be replaced by [f(x + ft) - f(x - h)]/(2h) with no change in the theory but with a dramatic improvement in the rate of convergence; that is, much more accurate estimates of /'(x) are obtained for a given value of ft. The kth derivative is the derivative of the (k - l)th derivative; they are de- k k (4) (5) noted by df/dx or /"(*), /'"(x), / (x), / (x),... Theorem 1.5 (Mean value theorem for derivatives) Let f(x) be continuous and differentiate in [a, 6] with \a\, |6| < oo. Then there is a point ζ in [a, b] so that ο — a Equivalently, given c and χ in [a, 6] there is a point ξ between χ and c so that f(x) = f(c) + f'(0(x-c) The special case of Theorem 1.5 with /(a) = f(b) = 0 is known as Rolle's theorem. It states that a if f() = fQ>) = 0> then there is a point ξ between a and b so that /'(£) = 0. This is derived from Theorem 1.5 by multiplying through by 6 - a, renaming a, b as x, c, and then applying the first form to the smaller interval [x, c] or [c, x], depending on the relation between χ and c. A very important tool in numerical analysis is the extension of the second part of Theorem 1.5 to use higher derivatives. Theorem 1.6 (Taylor's series with remainder) Let f(x) have η +1 continuous derivatives in [a, 6]. Given points χ and c in [a, b] we have Kx) = /(c) + /'(cX* - c) + / " ( O ^^ + / ' " ( c ) ^^ +... + / < " ) ( ) ^^ + C ß n + l) 0 r ζ! 3! ni where R \ has either one of the following forms (ξ is a point between χ and c): n+ (n + 1)! R i(x) = -^7 [\x-trfn+\lt)dt n+ If a function / depends on several variables, one can differentiate it with respect to one variable, say x, while keeping all the rest fixed. This is the partial derivative of / and is denoted by 6f/6x or f. Higher-order and mixed derivatives are defined by successive differentiation. Taylor's series x for functions of several variables is a direct extension of the formula in Theorem 1.6, although the number of terms in it grows rapidly. For two variables it is 2 /(*, V) = /(c, d) + f(x -c) + f(y-d)+ ^[f(x - cf + 2f(x - c)(y - d) + f (y - d)] + ... x y xx xy yy where all the partial derivatives are evaluated at the point (c, d). 4 MATHEMATICS AND COMPUTER SCIENCE BACKGROUND Theorem 1.7 (Chain rule for derivatives) Let f(x, y,..., z) have continuous first partial derivatives with respect to all its variables. Let χ = x(t), y = y(t),..., ζ = z(t) be continuously differentiate functions oft. Then g{t) = /(x(t), 2/(0, · · ·, z{t)) is continuously differentiable and g'(t) = fxx'(t) + fy'(t) + ... + fz'(t) y 2 Finally we state Theorem 1.8 (Fundamental theorem of algebra) Letpix) be a polynomial of degree η > 1, that is, 2 n p(x) = ao + a\x + a,2X +... + a x n where the ai are real or complex numbers and a ^ 0. Then there is a complex number ξ so that n p(O = 0. 1.1 PROBLEMS 2 2 1.1.1 Determine the following limits: (C) fix) - x = \ (D) x = cosß, fix) = sin θ _10_i (A) limî Zp^f 1.L5 Give the first four terms of the Taylor's series at c (B) lim ô χ—sîn χ for the following functions: x (C) \imx^/2 ^z^f (A) fix) = i\+x)k,c = 0 (D) lim _ n n+y/n (B) fix) = log 27, c = 1 (E) limô ψ- (C) /Or) = e*,c = 0 (F) lim _oo[log(n + 1) - log(n)] (D) f(x)= l/x,c= -1 n 1.1.2 Find the mean value ξ of Theorem 1.3 for the fol- (E) fix) = sin χ, c = π/6 lowing integrals: (F) fix) = \ / l - x\ c = -2 2 2 1.1.6 Obtain the Taylor's series for sin a; and cos χ with (A) jlxdx, with/Or) = s c = 0. Square each of them, add the results, and verify that 2 2 2 sin x + cos χ = 1. (B) Jq xdx, with/Or) = s 1.1.7 Find f and f for each of the following: x y 2 2 (C) J y/x~ix + \)dx with f(x) = Jx~ 2 2 0 (A) fix,y)=x -2xy+4y 2 2 2 (B) fix, y) = ix + y)/ix - y) (D) J^(x + \)dx with f(x) = x + 1 2 (C) fix,y)=x\0£iy/x) 1.1.3 Find the mean value ξ of Theorem 1.5 for the fol- (D) fix, y) = ln(y + z) where ^(x, 3/) = ye*'* lowing: 1.1.8 Find ci//cft for each of the following: 2 (A) fix) = x - Ax + 3, a = 1, b = 4 2 3 (A) fix, y) = x + xy, x = 3t, j/ = * (B)/(ar) = a?i, α = 0, 6 = 9 2 2 2 (B) fix,y)= y/x + y, x = t, y = e* - e~* iQ fix) = x + l/x, a = 1,6 = 9 1.1.9 Give the first six terms of the Taylor's series for the 1.1.4 Find /' and /" for each of the following: following: 3 (A) fix) = y/2x~+3 iA)fix,y) = ix + y), = 0,d = 0 c a;+1 (B) /(x) = sin2a: (B)/(x, /) = e /y, c = 0,d= 1 2 2 2 (C)/(^2/)=log0r + 2/), c=l,d = l

Numerical Methods in Software and Analysis

Similar Numerical Methods in Software and Analysis