f5-7

180 Chapter5. EvaluationofFunctionsSample 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).(the single real root when a, b, care real) and x2=−1 2(A+B)−a 3+i√ 3 2(A−B) x3=−1 2(A+B)−a 3−i√ 3 2(A−B)(5.6.18 ) (in that same case, a complex conjugate pair). Equations (5.6.13)–(5.6.16) are arrangedbothtominimizeroundofferror,andalso(aspointedoutbyA.J.Glassman)to ensure that no choice of branch for the complex cube root can result in the spurious loss of a distinct root. If you need to solve many cubic equations with only slightly different coeffi- cients, it is more efficient to use Newton’s method ( §9.4). CITED REFERENCES AND FURTHER READING: Weast,R.C.(ed.)1967, HandbookofTablesforMathematics , 3rded.(Cleveland:TheChemical Rubber Co.), pp. 130–133. Pachner, J. 1983, Handbookof NumericalAnalysis Applications (NewYork: McGraw-Hill), §6.1. McKelvey,J.P.1984, AmericanJournalofPhysics ,vol.52,pp.269–270;seealsovol.53,p.775, and vol. 55, pp. 374–375. 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/prime(x). Easy, right? The definition of the derivative, the limit as h→ 0of f/prime(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 fisfiercelyexpensive to compute, when you already have invested in computing f(x), and when, therefore, you want to get thederivativein nomorethan a singleadditionalfunctionevaluation. Insuch a situation, the remainingissue is to choose hproperly,an issue we nowdiscuss: Therearetwosourcesoferrorinequation(5.7.1),truncationerrorandroundoff error. Thetruncationerrorcomesfromhighertermsin theTaylorseries expansion, f(x+h)=f(x)+hf/prime(x)+1 2h2f/prime/prime(x)+1 6h3f/prime/prime/prime(x)+··· (5.7.2 ) whence f(x+h)−f(x) h=f/prime+1 2hf/prime/prime+··· (5.7.3 ) 5.7NumericalDerivatives 181Sample 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).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=1 0 .3and you blindly choose h=0.0001. Neither x=1 0 .3norx+h=1 0 .30001is a number with an exactrepresentationin binary;each is thereforerepresentedwith somefractional errorcharacteristicofthemachine’sfloating-pointformat, /epsilon1m,whosevalueinsingle precisionmaybe ∼10−7. Theerrorinthe effectivevalueof h,namelythedifference between x+handxasrepresentedinthemachine,isthereforeontheorderof /epsilon1mx, whichimpliesafractionalerrorin hoforder ∼/epsilon1mx/h∼10−2! Byequation(5.7.1) this immediatelyimplies at least the same largefractionalerrorin the derivative. WearriveatLesson1: Alwayschoose hsothat x+handxdifferbyanexactly representablenumber. This can usually be accomplishedby the programsteps temp =x+h h=temp−x(5.7.4 ) Some optimizing compilers, and some computers whose floating-point chips have higherinternalaccuracythanisstoredexternally,canfoilthistrick;ifso,itisusually enough to call a dummy subroutine donothing(temp) betweenthe two equations (5.7.4). This forces tempinto and out of addressable memory. With han “exact” number, the roundoff error in equation (5.7.1) is er∼ /epsilon1f|f(x)/h|. Here /epsilon1fis the fractional accuracy with which fis computed; for a simplefunctionthismaybecomparabletothemachineaccuracy, /epsilon1f≈/epsilon1m,butfora complicatedcalculation with additionalsources of inaccuracyit may be larger. The truncation error in equation (5.7.3) is on the order of et∼|hf/prime/prime(x)|. Varying hto minimize the sum er+etgives the optimal choice of h, h∼/radicalBigg /epsilon1ff f/prime/prime≈√/epsilon1fxc (5.7.5 ) where xc≡(f/f/prime/prime)1/2is the “curvature scale” of the function f, or “characteristic scale” over which it changes. In the absence of any other information, one often assumes xc=x(except near x=0where some other estimate of the typical x scale should be used). With the choice of equation (5.7.5), the fractional accuracy of the computed derivative is (er+et)/|f/prime|∼√/epsilon1f(ff/prime/prime/f/prime2)1/2∼√/epsilon1f (5.7.6 ) Here the last order-of-magnitude equality assumes that f,f/prime, and f/prime/primeall share the same characteristic length scale, usually the case. One sees that the simple finite-difference equation (5.7.1) gives at bestonly the square root of the machine accuracy /epsilon1m. If you can affordtwo functionevaluations for each derivativecalculation, then it is significantly better to use the symmetrized form f/prime(x)≈f(x+h)−f(x−h) 2h(5.7.7 ) 182 Chapter5. EvaluationofFunctionsSample 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).In this case, by equation (5.7.2), the truncation error is et∼h2f/prime/prime/prime. The roundoff error eris aboutthe same as before. The optimalchoiceof h, by a shortcalculation analogous to the one above, is now h∼/parenleftbigg/epsilon1ff f/prime/prime/prime/parenrightbigg1/3 ∼(/epsilon1f)1/3xc (5.7.8 ) and the fractional error is (er+et)/|f/prime|∼ (/epsilon1f)2/3f2/3(f/prime/prime/prime)1/3/f/prime∼(/epsilon1f)2/3(5.7.9 ) which will typically be an order of magnitude (single precision) or two orders of magnitude(doubleprecision) betterthanequation(5.7.6). WehavearrivedatLesson 2: Choose hto bethe correct powerof /epsilon1for/epsilon1mtimes a characteristicscale xc. You can easily derive the correct powers for other cases [1]. For a function of two dimensions, for example, and the mixed derivative formula ∂2f ∂x∂y=[f(x+h, y +h)−f(x+h, y−h)]−[f(x−h, y +h)−f(x−h, y−h)] 4h2 (5.7.10 ) the correct scaling is h∼/epsilon11/4 fxc. It is disappointing, certainly, that no simple finite-difference formula like equation(5.7.1)or(5.7.7)givesanaccuracycomparabletothemachineaccuracy /epsilon1m, oreventheloweraccuracytowhich fisevaluated, /epsilon1f. Aretherenobettermethods? Yes,thereare. All,however,involveexplorationofthefunction’sbehaviorover scalescomparableto xc,plussomeassumptionofsmoothness,oranalyticity,sothat thehigh-ordertermsinaTaylorexpansionlikeequation(5.7.2)havesomemeaning. Suchmethodsalso involvemultipleevaluationsofthe function f, so their increased accuracy must be weighed against increased cost. Thegeneralideaof“Richardson’sdeferredapproachtothelimit”isparticularly attractive. For numerical integrals, that idea leads to so-called Romberg integration(forreview, 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-differencecalculation to produce both an extrapolation of higher order, and also extrapolations of previous, lower,ordersbut with smaller scales h. Ridders [2]has givena nice implementation of this idea; the followingprogram, dfridr, is based on his algorithm,modified by animprovedterminationcriterion. Inputtotheroutineisafunction f(called func), a position x, and alargeststepsize h(more analogous to what we have called xc abovethanto whatwe havecalled h). Outputis thereturnedvalueofthederivative, and an estimate of its error, err. FUNCTION dfridr(func,x,h,err) INTEGER NTABREAL dfridr,err,h,x,func,CON,CON2,BIG,SAFEPARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.) EXTERNAL func C USES func Returns the derivative of a function func at a point xby Ridders’ method of polynomial extrapolation. The value his input as an estimated initial stepsize; it need not be small, 5.7NumericalDerivatives 183Sample 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).but rather should be an increment in xover which func changes substantially .A n e s t i m a t e of the error in the derivative is returned as err. Parameters: Stepsize is decreased by CON at each iteration. Max size of tableau is set by NTAB . Return when error is SAFE worse than the best so far. INTEGER i,j REAL errt,fac,hh,a(NTAB,NTAB) if(h.eq.0.) pause ’h must be nonzero in dfridr’hh=ha(1,1)=(func(x+hh)-func(x-hh))/(2.0*hh) err=BIG do 12i=2,NTAB Successive columns in the Neville tableau will go to smaller stepsizes and higher orders of extrapolation. hh=hh/CON a(1,i)=(func(x+hh)-func(x-hh))/(2.0*hh) Try new, smaller stepsize. fac=CON2 do11j=2,i Compute extrapolations of various orders, requiring no new function evaluations. a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.) fac=CON2*fac errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1))) The error strategy is to compare each new extrapolation to one order lower, both atthe present stepsize and the previous one. if (errt.le.err) then If error is decreased, save the improved answer. err=errtdfridr=a(j,i) endif enddo 11 if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err)return If higher order is worse by a significant factor SAFE , then quit early. enddo 12 return END Indfridr,thenumberofevaluationsof funcistypically6to12,butisallowed to be as great as 2 ×NTAB. As a function of input h, it is typical for the accuracy to getbetterashis made larger, until a sudden point is reached where nonsensical extrapolation produces early return with a large error. You should therefore choosea fairly large value for h, but monitor the returned value err, decreasing hif it is not small. For functions whose characteristic xscale is of order unity, we typically take hto be a few tenths. Besides Ridders’ method, there are other possible techniques. If your function is fairly smooth, and you know that you will want to evaluate its derivative manytimes at arbitrary points in some interval, then it makes sense to construct a Chebyshevpolynomialapproximationtothefunctioninthatinterval,andtoevaluate the derivative directly from the resulting Chebyshev coefficients. This method isdescribed in §§5.8–5.9, following. Another technique applies when the function consists of data that is tabulated at equally spaced intervals, and perhaps also noisy. One might then want, at each point, to least-squares fita polynomial of some degree M, using an additional number n Lof points to the left and some number nRof points to the right of each desired xvalue. The estimated derivative is then the derivative of the resulting fittedpolynomial. Averyefficientwaytodothisconstructionis viaSavitzky-Golay smoothing filters, which will be discussed later, in §14.8. There we will give a routineforgettingfiltercoefficientsthatnotonlyconstructthefittingpolynomialbut, in the accumulation of a single sum of data points times filter coefficients, evaluate it as well. In fact, the routine given, savgol, has an argument ldthat determines which derivative of the fitted polynomial is evaluated. For the first derivative, the 184 Chapter5. EvaluationofFunctionsSample 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).appropriate setting is ld=1, and the value of the derivative is the accumulated sum divided by the sampling interval h. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall), §§5.4–5.6. [1] Ridders, C.J.F. 1982, Advances in Engineering Software , vol. 4, no. 2, pp. 75–76. [2] 5.8 Chebyshev Approximation The Chebyshev polynomial of degree nis denoted Tn(x), and is given by the explicit formula Tn(x)=c o s ( narccos x)( 5.8.1 ) This may look trigonometric at first glance (and there is in fact a close relation between the Chebyshev polynomials and the discrete Fourier transform); however (5.8.1) can be combined with trigonometric identities to yield explicit expressionsforT n(x)(see Figure 5.8.1), T0(x)=1 T1(x)=x T2(x)=2 x2−1 T3(x)=4 x3−3x T4(x)=8 x4−8x2+1 ··· Tn+1(x)=2 xT n(x)−Tn−1(x)n≥1.(5.8.2 ) (There also exist inverse formulas for the powers of xin terms of the Tn’s — see equations 5.11.2-5.11.3.) TheChebyshevpolynomialsareorthogonalintheinterval [−1,1]overaweight (1−x2)−1/2. In particular, /integraldisplay1 −1Ti(x)Tj(x)√ 1−x2dx=/braceleftBigg0 i/negationslash=j π/2 i=j/negationslash=0 πi =j=0(5.8.3 ) The polynomial Tn(x)hasnzeros in the interval [−1,1], and they are located at the points x=c o s/parenleftbiggπ(k−1 2) n/parenrightbigg k=1,2,...,n (5.8.4 )