f9-1

9.1BracketingandBisection 343Sample 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).CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 5. Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapters 2, 7, and 14. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), Chapter 8. Householder, A.S. 1970, The Numerical Treatment of a Single Nonlinear Equation (New York: McGraw-Hill). 9.1 Bracketing and Bisection We will say that a root is bracketed in the interval (a, b )iff(a)andf(b)have opposite signs. If the function is continuous, then at least one root must lie in that interval (the intermediate value theorem ). If the function is discontinuous, but bounded, then instead of a root there might be a step discontinuity which crosses zero (see Figure 9.1.1). For numerical purposes, that might as well be a root, since the behavioris indistinguishablefromthe case of a continuousfunctionwhose zero crossing occurs in between two “adjacent” floating-point numbers in a machine’s finite-precision representation. Only for functions with singularities is there thepossibility that a bracketed root is not really there, as for example f(x)=1 x−c(9.1.1 ) Some root-finding algorithms (e.g., bisection in this section) will readily converge tocin (9.1.1). Luckily there is not much possibility of your mistaking c,o ra n y number xclose to it, for a root, since mere evaluation of |f(x)|will give a very large, rather than a very small, result. If you are given a function in a black box, there is no sure way of bracketing its roots, or of evendeterminingthat it has roots. If you like pathologicalexamples, thinkabouttheproblemoflocatingthetworealrootsofequation(3.0.1),whichdips below zero only in the ridiculously small interval of about x=π±10−667. In the next chapter we will deal with the related problem of bracketing a function’s minimum. There it is possible to give a procedure that always succeeds; in essence, “Go downhill, taking steps of increasing size, until your function startsbackuphill.” Thereisnoanalogousprocedureforroots. Theprocedure“godownhill until your function changes sign,” can be foiled by a function that has a simple extremum. Nevertheless, if you are prepared to deal with a “failure” outcome, this procedure is often a good first start; success is usual if your function has opposite signs in the limit x→± ∞. 344 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).a b (b)x1 ef c x 1dab ba(c) (d)(a) x2x3 Figure 9.1.1. Some situations encountered while root finding: (a) shows an isolated root x1bracketed by two points aand bat which the function has opposite signs; (b) illustrates that there is not necessarily a sign change in the function near a double root (in fact, there is not necessarily a root!); (c) is a pathological function with many roots; in (d) the function has opposite signs at points aand b, but the points bracket a singularity, not a root. 9.1BracketingandBisection 345Sample 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).SUBROUTINE zbrac(func,x1,x2,succes) INTEGER NTRY REAL x1,x2,func,FACTOR EXTERNAL funcPARAMETER (FACTOR=1.6,NTRY=50) Given a function func and an initial guessed range x1tox2, the routine expands the range geometrically until a root is bracketed by the returned values x1andx2(in which case succes returns as .true. ) or until the range becomes unacceptably large (in which case succes returns as .false. ). INTEGER j REAL f1,f2LOGICAL succesif(x1.eq.x2)pause ’you have to guess an initial range in zbrac’ f1=func(x1) f2=func(x2)succes=.true. do 11j=1,NTRY if(f1*f2.lt.0.)returnif(abs(f1).lt.abs(f2))then x1=x1+FACTOR*(x1-x2) f1=func(x1) else x2=x2+FACTOR*(x2-x1) f2=func(x2) endif enddo 11 succes=.false. return END Alternatively, you might want to “look inward ”on an initial interval, rather than“look outward ”from it, asking if there are any roots of the function f(x)in the interval from x1tox2when a search is carried out by subdivision into nequal intervals. The following subroutine returns brackets for up to nbdistinct intervals which each contain one or more roots. SUBROUTINE zbrak(fx,x1,x2,n,xb1,xb2,nb) INTEGER n,nbREAL x1,x2,xb1(nb),xb2(nb),fx EXTERNAL fx Given a function fxdefined on the interval from x1-x2 subdivide the interval into nequally spaced segments, and search for zero crossings of the function. nbis input as the maxi- mum number of roots sought, and is reset to the number of bracketing pairs xb1(1:nb) , xb2(1:nb) that are found. INTEGER i,nbbREAL dx,fc,fp,x nbb=0 x=x1dx=(x2-x1)/n Determine the spacing appropriate to the mesh. fp=fx(x) do 11i=1,n Loop over all intervals x=x+dxfc=fx(x) if(fc*fp.le.0.) then If a sign change occurs then record values for the bounds. nbb=nbb+1xb1(nbb)=x-dxxb2(nbb)=x if(nbb.eq.nb)goto 1 endiffp=fc enddo 11 346 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).1 continue nb=nbb return END Bisection Method Once we knowthat an interval contains a root, several classical proceduresare available to re fine it. These proceed with varying degrees of speed and sureness towardstheanswer. Unfortunately,themethodsthatareguaranteedtoconvergeplodalongmostslowly,whilethosethatrushtothesolutioninthebestcasescanalsodash rapidlytoin finitywithoutwarningifmeasuresarenottakentoavoidsuchbehavior. Thebisection method is one that cannot fail. It is thus not to be sneered at as a method for otherwise badly behaved problems. The idea is simple. Over some intervalthefunctionisknowntopassthroughzerobecauseitchangessign. Evaluatethe function at the interval ’s midpoint and examine its sign. Use the midpoint to replacewhicheverlimithasthesamesign. Aftereachiterationtheboundscontaining the root decrease by a factor of two. If after niterations the root is known to be within an interval of size /epsilon1 n, then after the next iteration it will be bracketed within an interval of size /epsilon1n+1=/epsilon1n/2( 9.1.2 ) neither more nor less. Thus, we know in advance the number of iterations required to achieve a given tolerance in the solution, n=l o g2/epsilon10 /epsilon1(9.1.3 ) where /epsilon10is the size of the initially bracketing interval, /epsilon1is the desired ending tolerance. Bisection mustsucceed. If the interval happens to contain two or more roots, bisectionwill findoneofthem. Iftheintervalcontainsnorootsandmerelystraddles a singularity, it will converge on the singularity. Whenamethodconvergesasafactor(lessthan1)timesthepreviousuncertainty tothefirstpower(asisthecaseforbisection),itissaidtoconverge linearly. Methods that converge as a higher power, /epsilon1n+1=constant ×(/epsilon1n)mm> 1( 9.1.4 ) are said to convergesuperlinearly. In other contexts “linear”convergencewould be termed“exponential, ”or“geometrical. ”Thatisnottoobadatall: Linearconvergence meansthatsuccessivesigni ficantfiguresarewonlinearlywith computationaleffort. It remains to discuss practical criteria for convergence. It is crucial to keep in mind that computers use a fixed number of binary digits to represent floating-point numbers. While your function might analytically pass through zero, it is possible that its computed value is never zero, for any floating-point argument. One must decide what accuracy on the root is attainable: Convergence to within 10−6in absolute value is reasonablewhen the root lies near 1, but certainly unachievableif 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