f9-2

9.2SecantMethod,FalsePositionMethod,andRidders’Method 347Sample 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 root lies near 1026. One might thus think to specify convergence by a relative (fractional) criterion, but this becomes unworkable for roots near zero. To be mostgeneral,theroutinesbelowwillrequireyoutospecifyanabsolutetolerance,suchthat iterations continue until the interval becomes smaller than this tolerance in absolute units. Usuallyyoumaywishtotakethetolerancetobe /epsilon1(|x 1|+|x2|)/2where /epsilon1isthe machineprecisionand x1andx2aretheinitialbrackets. Whentherootliesnearzero you ought to consider carefully what reasonable tolerance means for your function. The followingroutinequits after 40 bisections in anyevent, with 2−40≈10−12. FUNCTION rtbis(func,x1,x2,xacc) INTEGER JMAX REAL rtbis,x1,x2,xacc,funcEXTERNAL func PARAMETER (JMAX=40) Maximum allowed number of bisections. Using bisection, find the root of a function func known to lie between x1andx2.T h e root, returned as rtbis , will be refined until its accuracy is ±xacc . INTEGER j REAL dx,f,fmid,xmid fmid=func(x2)f=func(x1)if(f*fmid.ge.0.) pause ’root must be bracketed in rtbis’ if(f.lt.0.)then Orient the search so that f>0 lies at x+dx . rtbis=x1dx=x2-x1 else rtbis=x2dx=x1-x2 endif do 11j=1,JMAX Bisection loop. dx=dx*.5xmid=rtbis+dx fmid=func(xmid) if(fmid.le.0.)rtbis=xmidif(abs(dx).lt.xacc .or. fmid.eq.0.) return enddo 11 pause ’too many bisections in rtbis’END 9.2 Secant Method, False Position Method, and Ridders’ Method For functions that are smooth near a root, the methods known respectively asfalse position (orregula falsi ) andsecant method generally converge faster than bisection. In both of these methods the function is assumed to be approximately linearinthelocalregionofinterest,andthenextimprovementintherootistakenas the point where the approximatingline crosses the axis. After each iteration one ofthe previousboundarypointsis discardedin favorofthelatest estimateof theroot. Theonlydifference between the methods is that secant retains the most recent of the prior estimates (Figure 9.2.1; this requires an arbitrary choice on the first iteration),whilefalsepositionretainsthatpriorestimateforwhichthefunctionvalue 348 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).f(x)2 3 4 1x Figure 9.2.1. Secant method. Extrapolation or interpolation lines (dashed) are drawn through the two most recently evaluated points, whether or not they bracket the function. The points are numbered in the order that they are used. f(x) x432 1 Figure 9.2.2. False position method. Interpolation lines (dashed) are drawn through the most recent pointsthat bracket the root . In this example, point 1 thus remains “active”for many steps. False position converges less rapidly than the secant method, but it is more certain. 9.2SecantMethod,FalsePositionMethod,andRidders’Method 349Sample 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).2f(x) 134x Figure 9.2.3. Example where both the secant and false position methods will take many iterations to arrive at the true root. This function would be dif ficult for many other root- finding methods. has opposite sign from the function value at the current best estimate of the root, so that the two points continue to bracket the root (Figure 9.2.2). Mathematically, the secant method converges more rapidly near a root of a suf ficiently continuous function. Its order of convergencecan be shown to be the “golden ratio ”1.618..., so that lim k→∞|/epsilon1k+1|≈const×|/epsilon1k|1.618(9.2.1 ) Thesecant methodhas, however,the disadvantagethatthe rootdoesnot necessarily remain bracketed. For functions that are notsufficiently continuous, the algorithm can therefore not be guaranteed to converge: Local behavior might send it offtowards in finity. False position, since it sometimes keeps an older rather than newer function evaluation, has a lower order of convergence. Since the newer function value willsometimes be kept, the methodis often superlinear,but estimation of its exact order is not so easy. Here are sample implementations of these two related methods. While these methods are standard textbook fare, Ridders’ method , described below, or Brent’s method,inthenextsection,arealmostalwaysbetterchoices. Figure9.2.3showsthe behavior of secant and false-position methods in a dif ficult situation. FUNCTION rtflsp(func,x1,x2,xacc) INTEGER MAXIT REAL rtflsp,x1,x2,xacc,func EXTERNAL funcPARAMETER (MAXIT=30) Set to the maximum allowed number of iterations. 350 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).Using the false position method, find the root of a function func known to lie between x1 andx2. The root, returned as rtflsp , is refined until its accuracy is ±xacc . INTEGER j REAL del,dx,f,fh,fl,swap,xh,xlfl=func(x1) fh=func(x2) Be sure the interval brackets a root. if(fl*fh.gt.0.) pause ’root must be bracketed in rtflsp’if(fl.lt.0.)then Identify the limits so that xlcorresponds to the low side. xl=x1 xh=x2 else xl=x2xh=x1 swap=fl fl=fhfh=swap endif dx=xh-xldo 11j=1,MAXIT False position loop. rtflsp=xl+dx*fl/(fl-fh) Increment with respect to latest value. f=func(rtflsp) if(f.lt.0.) then Replace appropriate limit. del=xl-rtflsp xl=rtflsp fl=f else del=xh-rtflsp xh=rtflsp fh=f endif dx=xh-xl if(abs(del).lt.xacc.or.f.eq.0.)return Convergence. enddo 11 pause ’rtflsp exceed maximum iterations’ END FUNCTION rtsec(func,x1,x2,xacc) INTEGER MAXITREAL rtsec,x1,x2,xacc,funcEXTERNAL func PARAMETER (MAXIT=30) Maximum allowed number of iterations. Using the secant method, find the root of a function func thought to lie between x1and x2. The root, returned as rtsec , is refined until its accuracy is ±xacc . INTEGER j REAL dx,f,fl,swap,xl fl=func(x1)f=func(x2) if(abs(fl).lt.abs(f))then Pick the bound with the smaller function value as the most recent guess. rtsec=x1 xl=x2swap=fl fl=f f=swap else xl=x1 rtsec=x2 endifdo 11j=1,MAXIT Secant loop. dx=(xl-rtsec)*f/(f-fl) Increment with respect to latest value. xl=rtsecfl=f rtsec=rtsec+dx 9.2SecantMethod,FalsePositionMethod,andRidders’Method 351Sample 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).f=func(rtsec) if(abs(dx).lt.xacc.or.f.eq.0.)return Convergence. enddo 11 pause ’rtsec exceed maximum iterations’ END Ridders’Method A powerful variant on false position is due to Ridders [1]. When a root is bracketed between x1andx2, Ridders ’methodfirst evaluates the function at the midpoint x3=(x1+x2)/2. It then factors out that unique exponential function which turns the residual function into a straight line. Speci fically, it solves for a factor eQthat gives f(x1)−2f(x3)eQ+f(x2)e2Q=0 ( 9.2.2 ) This is a quadratic equation in eQ, which can be solved to give eQ=f(x3)+sign [f(x2)]/radicalbig f(x3)2−f(x1)f(x2) f(x2)(9.2.3 ) Now the false positionmethodis applied,not to the values f(x1),f(x3),f(x2),b u t to the values f(x1),f(x3)eQ,f(x2)e2Q, yielding a new guess for the root, x4. The overall updating formula (incorporating the solution 9.2.3) is x4=x3+(x3−x1)sign [f(x1)−f(x2)]f(x3)/radicalbig f(x3)2−f(x1)f(x2)(9.2.4 ) Equation (9.2.4) has some very nice properties. First, x4is guaranteed to lie in the interval (x1,x2), so the method never jumps out of its brackets. Second, the convergence of successive applications of equation (9.2.4) is quadratic , that is, m=2in equation (9.1.4). Since each application of (9.2.4) requires two function evaluations, the actual order of the method is√ 2, not 2; but this is still quite respectablysuperlinear: thenumberofsigni ficantdigitsintheanswerapproximately doubleswith each two functionevaluations. Third,takingout the function ’s“bend” via exponential (that is, ratio) factors, rather than via a polynomial technique (e.g., fitting a parabola), turns out to give an extraordinarily robust algorithm. In both reliabilityandspeed,Ridders ’methodisgenerallycompetitivewiththemorehighly developedandbetterestablished(butmorecomplicated)methodofVanWijngaarden, Dekker, and Brent, which we next discuss. FUNCTION zriddr(func,x1,x2,xacc) INTEGER MAXIT REAL zriddr,x1,x2,xacc,func,UNUSED PARAMETER (MAXIT=60,UNUSED=-1.11E30)EXTERNAL func C USES func Using Ridders’ method, return the root of a function func known to lie between x1and x2. The root, returned as zriddr , will be refined to an approximate accuracy xacc . INTEGER j REAL fh,fl,fm,fnew,s,xh,xl,xm,xnew 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.