f1-2

18 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).and routines psdesandran4in§7.5. We use a notation like Z’3F800000’ , which is consistent with the new FORTRAN-90 standard, but you may need to change this to, e.g., x’3f800000’ ,’3F800000’X ,o re v e n 16#3F800000 . Inextremis,youcan convertthe hexvaluesto decimalintegers;but notethat most compilerswill require anegativedecimalintegeras thevalueofa hexconstantwith its high-orderbitset. Asalreadymentionedin §1.0,thenotation a(1:m),inprogramcommentsandin thetext,denotesthearrayelementrange a(1),a(2),...,a(m). Likewise,notations likeb(2:7)orc(1:m,1:n) aretobeinterpretedasdenotingrangesofarrayindices. CITED REFERENCES AND FURTHER READING: Kernighan, B.W. 1978, The Elements of Programming Style (New York: McGraw-Hill). [1] Yourdon,E.1975, TechniquesofProgramStructureandDesign (EnglewoodCliffs,NJ:Prentice- Hall). [2] Meissner,L.P.andOrganick,E.I.1980, Fortran77FeaturingStructuredProgramming (Reading, MA: Addison-Wesley). [3] Hoare, C.A.R. 1981, Communications of the ACM , vol. 24, pp. 75–83. Wirth, N. 1983, Programming in Modula-2 , 3rd ed. (New York: Springer-Verlag). [4] Stroustrup, B. 1986, The C++Programming Language (Reading, MA: Addison-Wesley). [5] Borland International, Inc. 1989, Turbo Pascal 5.5 Object-Oriented Programming Guide (Scotts Valley, CA: Borland International). [6] Meeus, J. 1982, Astronomical Formulae for Calculators , 2nd ed., revised and enlarged (Rich- mond, VA: Willmann-Bell). [7] Hatcher,D.A.1984, QuarterlyJournaloftheRoyalAstronomicalSociety ,vol.25,pp.53–55;see alsoop. cit.1985, vol. 26, pp. 151–155, and 1986, vol. 27, pp. 506–507. [8] 1.2 Error, Accuracy, and Stability Althoughweassumenopriortrainingofthereaderinformalnumericalanalysis, we will need to presume a common understanding of a few key concepts. We will define these briefly in this section. Computers store numbers not with infinite precision but rather in some ap- proximation that can be packed into a fixed number of bits(binary digits) or bytes (groups of 8 bits). Almost all computers allow the programmer a choice among several different such representations ordata types . Data types can differ in the number of bits utilized (the wordlength ), but also in the more fundamental respect of whether the stored number is represented in fixed-point (also called integer)o r floating-point (also called real) format. A number in integer representation is exact. Arithmetic between numbers in integerrepresentationisalsoexact,withtheprovisosthat(i)theanswerisnotoutside the range of (usually, signed) integers that can be represented, and (ii) that division is interpretedas producinganintegerresult,throwingaway anyintegerremainder. 1.2Error,Accuracy,andStability 19Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).= =====1 100110 011000 011000 010000 011000 010000 110110 011000 000001 111010 101010 000000 001000 000000 001000 001000 000000 001000 000000 001000 001000 001000 001000 001000 001000 001000 000000 000000 001000 000000 001000 00000(a) (b) (c) (d) (e) (f) 1⁄2 3 1⁄4 10−7 3 + 10−7...23-bit mantissathis bit couldbe “phantom” sign bit 8-bit exponent Figure 1.2.1. Floating point representations of numbers in a typical 32-bit (4-byte) format. (a) The number 1/2(note the bias in the exponent); (b) the number 3; (c) the number 1/4; (d) the number 10−7, represented to machine accuracy; (e) the same number 10−7, but shifted so as to have the same exponent as the number 3; with this shifting, all signi ficance is lost and 10−7becomes zero; shifting to a common exponent must occur before two numbers can be added; (f) sum of the numbers 3+1 0−7, which equals 3tomachine accuracy. Eventhough 10−7can berepresented accurately byitself, itcannot accurately be added to a much larger number. Infloating-pointrepresentation,anumberisrepresentedinternallybyasignbit s(interpreted as plus or minus), an exact integer exponent e, and an exact positive integer mantissa M. Taken together these represent the number s×M×Be−E(1.2.1 ) where Bis the base of the representation (usually B=2, but sometimes B=1 6), andEis thebiasof the exponent, a fixed integer constant for any given machine and representation. An example is shown in Figure 1.2.1. Severalfloating-point bit patterns can represent the same number. If B=2, for example, a mantissa with leading (high-order)zero bits can be left-shifted, i.e., multipliedby a powerof 2, if the exponentis decreasedbya compensatingamount. Bit patterns that are “as left-shifted as they can be ”are termed normalized . Most computers always produce normalized results, since these don ’t waste any bits of the mantissa and thus allow a greater accuracy of the representation. Since the high-orderbitofaproperlynormalizedmantissa(when B=2)isalwaysone,some computers don ’t store this bit at all, giving one extra bit of signi ficance. Arithmeticamongnumbersin floating-pointrepresentationis notexact,evenif the operandshappento be exactlyrepresented(i.e., haveexactvalues in the formof equation1.2.1). For example, two floating numbersare added by first right-shifting (dividing by two) the mantissa of the smaller (in magnitude) one, simultaneously increasing its exponent,until the two operands have the same exponent. Low-order(least signi ficant) bits of the smaller operand are lost by this shifting. If the two operands differ too greatly in magnitude, then the smaller operand is effectively replaced by zero, since it is right-shifted to oblivion. The smallest (in magnitude) floating-point number which, when added to the floating-point number 1.0, produces a floating-point result different from 1.0 is termed the machine accuracy /epsilon1 m. A typical computer with B=2and a 32-bit wordlength has /epsilon1maround 3×10−8. (A more detailed discussion of machine 20 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).characteristics, and a program to determine them, is given in §20.1.) Roughly speaking,themachineaccuracy /epsilon1misthefractionalaccuracytowhich floating-point numbers are represented, corresponding to a change of one in the least signi ficant bit of the mantissa. Pretty much any arithmetic operation among floating numbers shouldbethoughtofasintroducinganadditionalfractionalerrorofatleast /epsilon1m. This type of error is called roundoff error . It is important to understand that /epsilon1mis not the smallest floating-point number thatcanberepresentedonamachine. Thatnumberdependsonhowmanybits there are inthe exponent,while /epsilon1mdependson howmanybits there arein the mantissa. Roundoff errors accumulate with increasing amounts of calculation. If, in the course of obtaining a calculated value, you perform Nsuch arithmetic operations, youmightbe so lucky as to have a total roundoff error on the order of√ N/epsilon1 m,i f the roundoff errors come in randomly up or down. (The square root comes from arandom-walk.) However,thisestimatecanbeverybadlyoffthemarkfortworeasons: (i) It very frequently happens that the regularities of your calculation, or the peculiaritiesofyourcomputer,causetheroundofferrorstoaccumulatepreferentiallyin one direction. In this case the total will be of order N/epsilon1 m. (ii) Some especially unfavorable occurrences can vastly increase the roundoff error of single operations. Generally these can be traced to the subtraction of two very nearly equal numbers, giving a result whose only signi ficant bits are those (few) low-order ones in which the operands differed. You might think that such a“coincidental ”subtraction is unlikely to occur. Not always so. Some mathematical expressionsmagnifyitsprobabilityofoccurrencetremendously. Forexample,inthe familiar formula for the solution of a quadratic equation, x=−b+√ b2−4ac 2a(1.2.2 ) the addition becomes delicate and roundoff-pronewhenever ac/lessmuchb2. (In§5.6 we will learn how to avoid the problem in this particular case.) Roundoff error is a characteristic of computer hardware. There is another, different, kind of error that is a characteristic of the program or algorithm used, independent of the hardware on which the program is executed. Many numericalalgorithms compute “discrete”approximations to some desired “continuous ”quan- tity. For example, an integral is evaluated numerically by computing a function at a discrete set of points, rather than at “every”point. Or, a function may be evaluated by summing a finite number of leading terms in its in finite series, rather than all in finity terms. In cases like this, there is an adjustable parameter, e.g., the number of points or of terms, such that the “true”answer is obtained only when that parameter goes to in finity. Any practical calculation is done with a finite, but sufficiently large, choice of that parameter. Thediscrepancybetweenthetrueanswerandtheanswerobtainedinapractical calculation is called the truncation error . Truncation error would persist even on a hypothetical, “perfect”computerthathadanin finitelyaccuraterepresentationandno roundofferror. As a general rule there is not much that a programmercan do about roundofferror, other than to choose algorithms that do not magnifyit unnecessarily (see discussionof “stability”below). Truncationerror,ontheotherhand,is entirely under the programmer ’s control. In fact, it is only a slight exaggeration to say 1.2Error,Accuracy,andStability 21Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).that clever minimization of truncation error is practically the entire content of the field of numerical analysis! Most of the time, truncation error and roundoff error do not strongly interact withoneanother. Acalculationcanbeimaginedashaving, first,thetruncationerror that it wouldhaveifrunonan in finite-precisioncomputer, “plus”the roundofferror associated with the number of operations performed. Sometimes, however, an otherwise attractive method can be unstable. This means that any roundofferror that becomes “mixed into ”the calculation at an early stageissuccessivelymagni fieduntilitcomestoswampthetrueanswer. Anunstable method would be useful on a hypothetical, perfect computer; but in this imperfectworld it is necessary for us to require that algorithms be stable —or if unstable that we use them with great caution. Here is a simple, if somewhat arti ficial, example of an unstable algorithm: Suppose that it is desired to calculate all integer powers of the so-called “Golden Mean,”the number given by φ≡√ 5−1 2≈0.61803398 ( 1.2.3 ) It turns out (you can easily verify) that the powers φnsatisfy a simple recursion relation, φn+1=φn−1−φn(1.2.4 ) Thus,knowingthe firsttwovalues φ0=1andφ1=0.61803398 ,wecansuccessively apply(1.2.4)performingonlyasinglesubtraction,ratherthanaslowermultiplicationbyφ, at each stage. Unfortunately,therecurrence(1.2.4)alsohas anothersolution,namelythevalue − 1 2(√ 5+1 ). Since the recurrence is linear, and since this undesired solution has magnitudegreaterthanunity,anysmalladmixtureofitintroducedbyroundofferrors will growexponentially. Ona typicalmachinewith 32-bitwordlength,(1.2.4)starts togivecompletelywronganswersbyabout n=1 6,atwhichpoint φnisdowntoonly 10−4. Therecurrence(1.2.4)is unstable,andcannotbeusedforthe purposestated. We will encounter the question of stability in many more sophisticated guises, later in this book. CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 1. Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 2. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §1.3. Wilkinson, J.H. 1964, Rounding Errors in Algebraic Processes (Englewood Cliffs, NJ: Prentice- Hall).