186 Chapter 5.Evaluation of Functions 5.7 Numerical Derivatives Imagine that you have a procedure which computes a function f(),and now you want to compute its derivative f'().Easy,right?The definition of the derivative,the limit as h-0 of f'()≈f+)-f) (5.7.1) h practically suggests the program:Pick a small value h;evaluate f(x+h);you probably have f()already evaluated,but if not,do it too;finally apply equation (5.7.1).What more needs to be said? 鱼 Quite a lot,actually.Applied uncritically,the above procedure is almost 二 guaranteed to produce inaccurate results.Applied properly,it can be the right way to compute a derivative only when the function f is fiercely expensive to compute, when you already have invested in computing f(z),and when,therefore,you want RECIPES to get the derivative in no more than a single additional function evaluation.In such a situation,the remaining issue is to choose h properly,an issue we now discuss: 言O。 9 There are two sources of error in equation(5.7.1).truncation error and roundoff error.The truncation error comes from higher terms in the Taylor series expansion, fe+)=f回+hf回+fa回)+f”a)+ (5.7.2) 。2 whence fe+)-f@=f+hf”+… 1 h (5.7.3) The roundoff error has various contributions.First there is roundoff error in h: Suppose,by way of an example,that you are at a point x=10.3 and you blindly choose h =0.0001.Neither x 10.3 nor x+h =10.30001 is a number with an exact representation in binary;each is therefore represented with some fractional error characteristic of the machine's floating-point format,em,whose value in single Recipes 10621 precision may be~10-7.The error in the effective value of h,namely the difference Numerica between x+h and x as represented in the machine,is therefore on the order of em, 4310 which implies a fractional error in h oforder~em/h10-2!By equation (5.7.1) Recipes this immediately implies at least the same large fractional error in the derivative. We arrive at Lesson 1:Always choose h so that +h and z differ by an exactly representable number.This can usually be accomplished by the program steps North temp=x+h (5.7.4) h temp- Some optimizing compilers,and some computers whose floating-point chips have higher internal accuracy than is stored externally,can foil this trick;if so,it is usually enough to declare temp as volatile,or else to call a dummy function donothing(temp)between the two equations(5.7.4).This forces temp into and out of addressable memory186 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 5.7 Numerical Derivatives Imagine that you have a procedure which computes a function f(x), and now you want to compute its derivative f (x). Easy, right? The definition of the derivative, the limit as h → 0 of f (x) ≈ f(x + h) − f(x) h (5.7.1) practically suggests the program: Pick a small value h; evaluate f(x + h); you probably have f(x) already evaluated, but if not, do it too; finally apply equation (5.7.1). What more needs to be said? Quite a lot, actually. Applied uncritically, the above procedure is almost guaranteed to produce inaccurate results. Applied properly, it can be the right way to compute a derivative only when the function f is fiercely expensive to compute, when you already have invested in computing f(x), and when, therefore, you want to get the derivative in no more than a single additional function evaluation. In such a situation, the remaining issue is to choose h properly, an issue we now discuss: There are two sources of error in equation (5.7.1), truncation error and roundoff error. The truncation error comes from higher terms in the Taylor series expansion, f(x + h) = f(x) + hf (x) + 1 2 h2f(x) + 1 6 h3f(x) + ··· (5.7.2) whence f(x + h) − f(x) h = f + 1 2 hf + ··· (5.7.3) The roundoff error has various contributions. First there is roundoff error in h: Suppose, by way of an example, that you are at a point x = 10.3 and you blindly choose h = 0.0001. Neither x = 10.3 nor x + h = 10.30001 is a number with an exact representation in binary; each is therefore represented with some fractional error characteristic of the machine’s floating-point format,  m, whose value in single precision may be ∼ 10−7. The error in the effective value of h, namely the difference between x + h and x as represented in the machine, is therefore on the order of  mx, which implies a fractional error in h of order ∼ mx/h ∼ 10−2! By equation (5.7.1) this immediately implies at least the same large fractional error in the derivative. We arrive at Lesson 1: Always choose h so that x+ h and x differ by an exactly representable number. This can usually be accomplished by the program steps temp = x + h h = temp − x (5.7.4) Some optimizing compilers, and some computers whose floating-point chips have higher internal accuracy than is stored externally, can foil this trick; if so, it is usually enough to declare temp as volatile, or else to call a dummy function donothing(temp) between the two equations (5.7.4). This forces temp into and out of addressable memory.