f9-3

352 Chapter9. RootFindingandNonlinearSetsof EquationsSample 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).fl=func(x1) fh=func(x2) if((fl.gt.0..and.fh.lt.0.).or.(fl.lt.0..and.fh.gt.0.))then xl=x1xh=x2 zriddr=UNUSED Any highly unlikely value, to simplify logic below. do 11j=1,MAXIT xm=0.5*(xl+xh)fm=func(xm) First of two function evaluations per it- eration. s=sqrt(fm**2-fl*fh) if(s.eq.0.)returnxnew=xm+(xm-xl)*(sign(1.,fl-fh)*fm/s) Updating formula. if (abs(xnew-zriddr).le.xacc) return zriddr=xnew fnew=func(zriddr) Second of two function evaluations per iteration. if (fnew.eq.0.) return if(sign(fm,fnew).ne.fm) then Bookkeeping to keep the root bracketed on next iteration. xl=xm fl=fmxh=zriddr fh=fnew else if(sign(fl,fnew).ne.fl) then xh=zriddr fh=fnew else if(sign(fh,fnew).ne.fh) then xl=zriddrfl=fnew else pause ’never get here in zriddr’ endif if(abs(xh-xl).le.xacc) return enddo 11 pause ’zriddr exceed maximum iterations’ else if (fl.eq.0.) then zriddr=x1 else if (fh.eq.0.) then zriddr=x2 else pause ’root must be bracketed in zriddr’ endifreturn END CITED REFERENCES AND FURTHER READING: Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill),§8.3. Ostrowski, A.M. 1966, Solutions of Equations and Systems of Equations , 2nd ed. (New York: Academic Press), Chapter 12. Ridders,C.J.F.1979, IEEETransactionsonCircuitsandSystems , vol.CAS-26,pp.979–980.[1] 9.3 Van Wijngaarden–Dekker–Brent Method While secant and false position formally converge faster than bisection, one finds in practice pathologicalfunctions for which bisection convergesmore rapidly. 9.3VanWijngaarden–Dekker–BrentMethod 353Sample 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).These can be choppy, discontinuous functions, or even smooth functions if the secondderivativechangessharplyneartheroot. Bisectionalwayshalvestheinterval,while secant and false position can sometimes spend many cycles slowly pulling distant bounds closer to a root. Ridders’ method does a much better job, but it too can sometimes be fooled. Is there a way to combine superlinear convergencewith the sureness of bisection? Yes. We cankeeptrackofwhethera supposedlysuperlinearmethodis actually converging the way it is supposed to, and, if it is not, we can intersperse bisection steps so as to guarantee at leastlinear convergence. This kind of super-strategy requires attention to bookkeeping detail, and also careful consideration of howroundofferrors can affect the guiding strategy. Also, we must be able to determine reliably when convergence has been achieved. Anexcellentalgorithmthatpayscloseattentiontothesematterswasdeveloped in the 1960s by van Wijngaarden, Dekker, and others at the Mathematical Center in Amsterdam, and later improved by Brent [1]. For brevity, we refer to the final form of the algorithm as Brent’s method . The method is guaranteed (by Brent) to converge, so long as the function can be evaluated within the initial interval known to contain a root. Brent’s method combines root bracketing, bisection, and inverse quadratic interpolation toconvergefromthe neighborhoodofa zerocrossing. While the false position and secant methods assume approximately linear behavior between twoprior root estimates, inverse quadratic interpolation uses three prior points to fit an inverse quadratic function ( xas a quadratic function of y) whose value at y=0is takenas thenextestimate of the root x. Of courseone musthavecontingencyplans for what to do if the root falls outside of the brackets. Brent’s method takes care of all that. Ifthethreepointpairsare [a, f (a)],[b, f (b)],[c, f (c)]thentheinterpolation formula (cf. equation 3.1.1) is x=[y−f(a)][y−f(b)]c [f(c)−f(a)][f(c)−f(b)]+[y−f(b)][y−f(c)]a [f(a)−f(b)][f(a)−f(c)] +[y−f(c)][y−f(a)]b [f(b)−f(c)][f(b)−f(a)](9.3.1 ) Setting yto zero givesa result forthe nextroot estimate, which canbe written as x=b+P/Q (9.3.2 ) where, in terms of R≡f(b)/f(c),S ≡f(b)/f(a),T ≡f(a)/f(c)( 9.3.3 ) we have P=S[T(R−T)(c−b)−(1−R)(b−a)] ( 9.3.4 ) Q=(T−1)(R−1)(S−1) ( 9.3.5 ) In practice bis the current best estimate of the root and P/Qought to be a “small” correction. Quadraticmethodsworkwellonlywhenthefunctionbehavessmoothly; 354 Chapter9. RootFindingandNonlinearSetsofEquationsSample 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).they run the serious risk of giving very bad estimates of the next root or causing machine failure by an inappropriate division by a very small number ( Q≈0). Brent’s method guards against this problem by maintaining brackets on the root and checking where the interpolation would land before carrying out the division. When the correction P/Qwould not land within the bounds, or when the bounds are not collapsing rapidly enough, the algorithm takes a bisection step. Thus, Brent’s methodcombines the sureness of bisection with the speed of a higher-order method when appropriate. We recommend it as the method of choice for general one-dimensional root finding where a function’s values only (and not its derivative or functional form) are available. FUNCTION zbrent(func,x1,x2,tol) INTEGER ITMAXREAL zbrent,tol,x1,x2,func,EPS EXTERNAL func PARAMETER (ITMAX=100,EPS=3.e-8) Using Brent’s method, find the root of a function func known to lie between x1andx2. The root, returned as zbrent , will be refined until its accuracy is tol. Parameters: Maximum allowed number of iterations, and machine floating-point precision. INTEGER iterREAL a,b,c,d,e,fa,fb,fc,p,q,r, * s,tol1,xm a=x1b=x2fa=func(a) fb=func(b) if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) * pause ’root must be bracketed for zbrent’ c=b fc=fbdo 11iter=1,ITMAX if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then c=a Rename a,b,cand adjust bounding interval d. fc=fad=b-a e=d endifif(abs(fc).lt.abs(fb)) then a=b b=c c=afa=fbfb=fc fc=fa endiftol1=2.*EPS*abs(b)+0.5*tol Convergence check. xm=.5*(c-b) if(abs(xm).le.tol1 .or. fb.eq.0.)then zbrent=breturn endif if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then s=fb/fa Attempt inverse quadratic interpolation. if(a.eq.c) then p=2.*xm*sq=1.-s else q=fa/fc r=fb/fcp=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) 9.4Newton-RaphsonMethodUsingDerivative 355Sample 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).endif if(p.gt.0.) q=-q Check whether in bounds. p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d Accept interpolation. d=p/q else d=xm Interpolation failed, use bisection. e=d endif else Bounds decreasing too slowly, use bisection. d=xme=d endif a=b Move last best guess to a. fa=fb if(abs(d) .gt. tol1) then Evaluate new trial root. b=b+d else b=b+sign(tol1,xm) endif fb=func(b) enddo 11 pause ’zbrent exceeding maximum iterations’zbrent=breturnEND CITED REFERENCES AND FURTHER READING: Brent,R.P.1973, AlgorithmsforMinimizationwithoutDerivatives (EnglewoodCliffs,NJ:Prentice- Hall), Chapters 3, 4. [1] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §7.2. 9.4 Newton-Raphson Method Using Derivative Perhapsthemostcelebratedofallone-dimensionalroot-findingroutinesis New- ton’smethod ,alsocalledthe Newton-Raphsonmethod . Thismethodisdistinguished from the methods of previous sections by the fact that it requires the evaluation of both the function f(x),andthe derivative f/prime(x), at arbitrary points x. The Newton-Raphson formula consists geometrically of extending the tangent line at a currentpoint xiuntil it crosses zero,thensetting thenextguess xi+1to theabscissa of that zero-crossing(see Figure 9.4.1). Algebraically, the method derives from thefamiliar Taylor series expansionof a functionin the neighborhoodof a point, f(x+δ)≈f(x)+f /prime(x)δ+f/prime/prime(x) 2δ2+.... (9.4.1 ) For small enough values of δ, and for well-behaved functions, the terms beyond linear are unimportant, hence f(x+δ)=0implies δ=−f(x) f/prime(x). (9.4.2 )