f9-0

Sample 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).Chapter 9. Root Finding and Nonlinear SetsofEquations 9.0 Introduction Wenowconsiderthatmostbasicoftasks,solvingequationsnumerically. While most equations are born with both a right-hand side and a left-hand side, one traditionally moves all terms to the left, leaving f(x)=0 ( 9.0.1 ) whosesolutionorsolutionsaredesired. Whenthereisonlyoneindependentvariable, the problemis one-dimensional ,namelyto find the root or roots of a function. With more than one independent variable, more than one equation can be satisfied simultaneously. You likely once learned the implicit function theorem which (in this context) gives us the hope of satisfying Nequations in Nunknowns simultaneously. Note that we have only hope, not certainty. A nonlinear set of equations may have no (real) solutions at all. Contrariwise, it may have more thanone solution. The implicit function theorem tells us that “generically”the solutions will be distinct, pointlike, and separated from each other. If, however, life is so unkind as to present you with a nongeneric, i.e., degenerate, case, then you can geta continuous family of solutions. In vector notation, we want to find one or more N-dimensional solution vectors xsuch that f(x)=0 (9.0.2 ) wherefis the N-dimensional vector-valued function whose components are the individual equations to be satisfied simultaneously. Don’t be fooled by the apparent notational similarity of equations (9.0.2) and (9.0.1). Simultaneous solutionof equationsin Ndimensionsis muchmore difficult thanfindingrootsintheone-dimensionalcase. Theprincipaldifferencebetweenoneandmanydimensionsisthat,inonedimension,itispossibletobracketor“trap”aroot between bracketingvalues, and then huntit down like a rabbit. In multidimensions, you can never be sure that the root is there at all until you have found it. Except in linear problems, root finding invariably proceeds by iteration, and this is equally true in one or in many dimensions. Starting from some approximatetrialsolution,ausefulalgorithmwill improvethesolutionuntilsomepredetermined convergencecriterion is satisfied. For smoothly varyingfunctions, good algorithms 340 9.0Introduction 341Sample 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).willalwaysconverge, providedthattheinitialguessis goodenough. Indeedonecan even determine in advance the rate of convergenceof most algorithms. Itcannotbeoveremphasized,however,howcruciallysuccessdependsonhaving a good first guess for the solution, especially for multidimensional problems. This crucialbeginningusuallydependsonanalysisratherthannumerics. Carefullycraftedinitial estimates reward you not only with reduced computational effort, but also with understanding and increased self-esteem. Hamming’s motto, “the purpose of computing is insight, not numbers,” is particularly apt in the area of finding roots. Youshouldrepeatthismottoaloudwheneveryourprogramconverges,withten-digit accuracy, to the wrong root of a problem, or whenever it fails to converge becausethere is actually noroot, or because there is a root but your initial estimate was not sufficiently close to it. “This talk of insight is all very well, but what do I actually do?” For one- dimensional root finding, it is possible to give some straightforward answers: You should try to get some idea of what your function looks like before trying to find its roots. If you need to mass-produce roots for many different functions, then youshould at least know what some typical members of the ensemble look like. Next, you should always bracket a root, that is, know that the function changes sign in an identified interval, before trying to converge to the root’s value. Finally (this is advice with which some daring souls might disagree, but we giveitnonetheless)neverletyouriterationmethodgetoutsideofthebestbracketingboundsobtainedat anystage. We will see belowthat somepedagogicallyimportant algorithms,suchas secantmethod orNewton-Raphson ,canviolatethislastconstraint, and are thus not recommended unless certain fixups are implemented. Multiple roots, or very close roots, are a real problem, especially if the multiplicity is an even number. In that case, there may be no readily apparent sign change in the function, so the notion of bracketing a root — and maintaining the bracket — becomes difficult. We are hard-liners: we nevertheless insist on bracketing a root, even if it takes the minimum-searchingtechniques of Chapter 10to determine whether a tantalizing dip in the function really does cross zero or not. (You can easily modify the simple golden section routine of §10.1 to return early if it detects a sign change in the function. And, if the minimum of the function isexactly zero, then you have found a doubleroot.) Asusual,wewanttodiscourageyoufromusingroutinesasblackboxeswithout understanding them. However, as a guide to beginners, here are some reasonable starting points: •Brent’s algorithm in §9.3 is the method of choice to find a bracketed root of a general one-dimensional function, when you cannot easily compute the function’s derivative. Ridders’ method ( §9.2) is concise, and a close competitor. •When you can compute the function’s derivative, the routine rtsafein §9.4,whichcombinesthe Newton-Raphsonmethodwith somebookkeep- ingonbounds,is recommended. Again,youmustfirst bracketyourroot. •Roots of polynomials are a special case. Laguerre’s method, in §9.5, is recommended as a starting point. Beware: Some polynomials areill-conditioned! •Finally, for multidimensional problems, the only elementary method is Newton-Raphson ( §9.6), which works verywell if you can supply a 342 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).good first guess of the solution. Try it. Then read the more advanced materialin §9.7forsomemorecomplicated,butgloballymoreconvergent, alternatives. Avoiding implementations for specific computers, this book must generally steer clear of interactive or graphics-related routines. We make an exception rightnow. Thefollowingroutine,whichproducesa crudefunctionplotwith interactively scaled axes, can save you a lot of grief as you enter the world of root finding. SUBROUTINE scrsho(fx) INTEGER ISCR,JSCR REAL fx EXTERNAL fxPARAMETER (ISCR=60,JSCR=21) Number of horizontal and vertical positions in display. For interactive CRT terminal use. Produce a crude graph of the function fxover the prompted-for interval x1,x2 . Query for another plot until the user signals satisfaction. INTEGER i,j,jzREAL dx,dyj,x,x1,x2,ybig,ysml,y(ISCR)CHARACTER*1 scr(ISCR,JSCR),blank,zero,yy,xx,ff SAVE blank,zero,yy,xx,ff DATA blank,zero,yy,xx,ff/’ ’,’-’,’l’,’-’,’x’/ 1 continue write (*,*) ’ Enter x1,x2 (= to stop)’ Query for another plot, quit if x1=x2 . read (*,*) x1,x2if(x1.eq.x2) returndo 11j=1,JSCR Fill vertical sides with character ’ l’. scr(1,j)=yy scr(ISCR,j)=yy enddo 11 do13i=2,ISCR-1 scr(i,1)=xx Fill top, bottom with character ’ -’. scr(i,JSCR)=xxdo 12j=2,JSCR-1 Fill interior with blanks. scr(i,j)=blank enddo 12 enddo 13 dx=(x2-x1)/(ISCR-1)x=x1ybig=0. Limits will include 0. ysml=ybig do 14i=1,ISCR Evaluate the function at equal intervals. Find the largest and smallest values. y(i)=fx(x) if(y(i).lt.ysml) ysml=y(i)if(y(i).gt.ybig) ybig=y(i) x=x+dx enddo 14 if(ybig.eq.ysml) ybig=ysml+1. Be sure to separate top and bottom. dyj=(JSCR-1)/(ybig-ysml) jz=1-ysml*dyj Note which row corresponds to 0. do15i=1,ISCR Place an indicator at function height and 0. scr(i,jz)=zero j=1+(y(i)-ysml)*dyj scr(i,j)=ff enddo 15 write (*,’(1x,1pe10.3,1x,80a1)’) ybig,(scr(i,JSCR),i=1,ISCR)do 16j=JSCR-1,2,-1 Display. write (*,’(12x,80a1)’) (scr(i,j),i=1,ISCR) enddo 16 write (*,’(1x,1pe10.3,1x,80a1)’) ysml,(scr(i,1),i=1,ISCR)write (*,’(12x,1pe10.3,40x,e10.3)’) x1,x2goto 1 END 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→± ∞.