188 Chapter 5. Evaluation of Functions the correct scaling ish It is disappointing,certainly,that no simple finite-difference formula like equation (5.7.1)or(5.7.7)gives an accuracy comparable to the machine accuracy em,or even the lower accuracy to which f is evaluated,ef.Are there no better methods? Yes.there are.All.however,involve exploration of the function's behavior over scales comparable to e,plus some assumption of smoothness,or analyticity,so that the high-order terms in a Taylor expansion like equation(5.7.2)have some meaning. Such methods also involve multiple evaluations of the function f.so their increased accuracy must be weighed against increased cost. 8 The general idea of"Richardson's deferred approach to the limit"is particularly attractive.For numerical integrals,that idea leads to so-called Romberg integration (for review.see 84.3).For derivatives,one seeks to extrapolate,to h-0,the result of finite-difference calculations with smaller and smaller finite values of h.By the use of Neville's algorithm(83.1),one uses each new finite-difference calculation to produce both an extrapolation of higher order,and also extrapolations of previous, ≌cM lower,orders but with smaller scales h.Ridders [2]has given a nice implementation of this idea;the following program,dfridr,is based on his algorithm,modified by an improved termination criterion.Input to the routine is a function f(called func), a position z,and a largest stepsize h(more analogous to what we have called above than to what we have called h).Output is the returned value of the derivative, 公 Press. and an estimate of its error,err. 9 #include <math.h> SCIENTIFIC #include "nrutil.h" #define CON 1.4 Stepsize is decreased by CON at each iteration. #define CON2 (CON*CON) #define BIG 1.0e30 #define NTAB 10 Sets maximum size of tableau. #define SAFE 2.0 Return when error is SAFE worse than the best so far. COMPUTING (ISBN 18881920 float dfridr(float (*func)(float),float x,float h,float terr) Returns the derivative of a function func at a point x by Ridders'method of polynomial extrapolation.The value h is input as an estimated initial stepsize;it need not be small,but 10621 rather should be an increment in x over which func changes substanrially.An estimate of the error in the derivative is returned as err. idge.org Numerical Recipes 43108 int i,ji float errt,fac,hh,**a,ans; (outside 膜 if (h ==0.0)nrerror("h must be nonzero in dfridr."); a=matrix(1,NTAB,1,NTAB); oftware. hh=hi a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); *err=BIG; for(i=2;1<=NTAB;1++)[ Successive columns in the Neville tableau will go to smaller stepsizes and higher orders of extrapolation. hh /CON; a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); Try new,smaller step fac=CON2; size. for(j=2;j<=1;j+){ Compute extrapolations of various orders,requiring a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0); no new function eval- fac=CON2*fac: uations. errt=FMAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));188 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). the correct scaling is h ∼  1/4 f xc. It is disappointing, certainly, that no simple finite-difference formula like equation (5.7.1) or (5.7.7) gives an accuracy comparable to the machine accuracy m, or even the lower accuracy to which f is evaluated,  f . Are there no better methods? Yes, there are. All, however, involve exploration of the function’s behavior over scales comparable to xc, plus some assumption of smoothness, or analyticity, so that the high-order terms in a Taylor expansion like equation (5.7.2) have some meaning. Such methods also involve multiple evaluations of the function f, so their increased accuracy must be weighed against increased cost. The general idea of “Richardson’s deferred approach to the limit” is particularly attractive. For numerical integrals, that idea leads to so-called Romberg integration (for review, see §4.3). For derivatives, one seeks to extrapolate, to h → 0, the result of finite-difference calculations with smaller and smaller finite values of h. By the use of Neville’s algorithm (§3.1), one uses each new finite-difference calculation to produce both an extrapolation of higher order, and also extrapolations of previous, lower, orders but with smaller scales h. Ridders [2] has given a nice implementation of this idea; the following program, dfridr, is based on his algorithm, modified by an improved termination criterion. Input to the routine is a function f (called func), a position x, and a largest stepsize h (more analogous to what we have called xc above than to what we have called h). Output is the returned value of the derivative, and an estimate of its error, err. #include <math.h> #include "nrutil.h" #define CON 1.4 Stepsize is decreased by CON at each iteration. #define CON2 (CON*CON) #define BIG 1.0e30 #define NTAB 10 Sets maximum size of tableau. #define SAFE 2.0 Return when error is SAFE worse than the best so far. float dfridr(float (*func)(float), float x, float h, float *err) Returns the derivative of a function func at a point x by Ridders’ method of polynomial extrapolation. The value h is input as an estimated initial stepsize; it need not be small, but rather should be an increment in x over which func changes substantially. An estimate of the error in the derivative is returned as err. { int i,j; float errt,fac,hh,**a,ans; if (h == 0.0) nrerror("h must be nonzero in dfridr."); a=matrix(1,NTAB,1,NTAB); hh=h; a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); *err=BIG; for (i=2;i<=NTAB;i++) { Successive columns in the Neville tableau will go to smaller stepsizes and higher orders of extrapolation. hh /= CON; a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); Try new, smaller step￾fac=CON2; size. for (j=2;j<=i;j++) { Compute extrapolations of various orders, requiring no new function eval￾uations. a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0); fac=CON2*fac; errt=FMAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));