f9-6

372 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).This equation, if used with iranging over the roots already polished, will prevent a tentative root from spuriously hopping to another one’s true root. It is an exampleof so-called zero suppression as an alternative to true deflation. Muller’smethod,whichwasdescribedabove,canalsobeusefulatthepolishing stage. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 7. [1] PetersG.,andWilkinson,J.H.1971, JournaloftheInstituteofMathematicsanditsApplications , vol. 8, pp. 16–35. [2] IMSL Math/Library UsersManual (IMSL Inc., 2500CityWest Boulevard, HoustonTX77042).[3] Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §8.9–8.13. [4] Adams, D.A. 1967, Communications of the ACM , vol. 10, pp. 655–658. [5] Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §4.4.3. [6] Henrici, P. 1974, Applied and Computational Complex Analysis , vol. 1 (New York: Wiley). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §§5.5–5.9. 9.6 Newton-Raphson Method for Nonlinear Systems of Equations Wemakeanextreme,butwhollydefensible,statement: Thereare nogood,gen- eralmethodsforsolvingsystemsofmorethanonenonlinearequation. Furthermore, itis nothardtosee why(verylikely)there neverwill be anygood,generalmethods: Consider the case of two dimensions, where we want to solve simultaneously f(x, y )=0 g(x, y )=0(9.6.1 ) The functions fandgare two arbitrary functions, each of which has zero contourlinesthatdividethe (x, y )planeintoregionswheretheirrespectivefunction is positive or negative. These zero contour boundaries are of interest to us. The solutionsthatweseekarethosepoints(ifany)thatarecommontothezerocontours offandg(see Figure9.6.1). Unfortunately,the functions fandghave,in general, norelationtoeachotheratall! Thereis nothingspecialabouta commonpointfromeither f’s point of view, or from g’s. In order to find all common points, which are thesolutionsofournonlinearequations,wewill(ingeneral)havetodoneithermore nor less than map out the full zero contours of both functions. Note further thatthe zero contours will (in general) consist of an unknownnumberof disjoint closed curves. Howcanweeverhopetoknowwhenwehavefoundallsuchdisjointpieces? For problems in more than two dimensions, we need to find points mutually commonto Nunrelatedzero-contourhypersurfaces,eachofdimension N−1.Y o u 9.6Newton-RaphsonMethodforNonlinearSystems ofEquations 373Sample 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).g=0 g=0f=0f=0f pos Mg pos f pos f pos f neg g=0g negg pos g neg g posy xno root here! two roots here Figure 9.6.1. Solution of two nonlinear equations in two unknowns. Solid curves refer to f(x, y ), dashed curves to g(x, y ). Each equation divides the (x, y )plane into positive and negative regions, bounded by zero curves. The desired solutions are the intersections of these unrelated zero curves. The number of solutions is a prioriunknown. see that root finding becomes virtually impossible without insight! You will almost always have to use additional information, speci fic to your particular problem, to answersuchbasicquestionsas, “DoIexpectauniquesolution? ”and“Approximately where?”Acton[1]has a good discussion of some of the particular strategies that can be tried. In this section we will discuss the simplest multidimensional root finding method, Newton-Raphson. This method gives you a very ef ficient means of converging to a root, if you have a suf ficiently good initial guess. It can also spectacularly fail to converge, indicating (though not proving) that your putativeroot does not exist nearby. In §9.7 we discuss more sophisticated implementations of the Newton-Raphson method, which try to improve on Newton-Raphson ’s poor globalconvergence. Amultidimensionalgeneralizationofthesecantmethod,calledBroyden’s method, is also discussed in §9.7. Atypicalproblemgives Nfunctionalrelationstobezeroed,involvingvariables x i,i=1,2,...,N: Fi(x1,x2,...,x N)=0 i=1,2,...,N. (9.6.2 ) We letxdenote the entire vector of values xiandFdenote the entire vector of functions Fi. In the neighborhood of x, each of the functions Fican be expanded in Taylor series Fi(x+δx)=Fi(x)+N/summationdisplay j=1∂F i ∂x jδx j+O(δx2). (9.6.3 ) 374 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).Thematrixofpartialderivativesappearinginequation(9.6.3)isthe Jacobian matrixJ: Jij≡∂F i ∂x j. (9.6.4 ) In matrix notation equation (9.6.3) is F(x+δx)=F(x)+J·δx+O(δx2). (9.6.5 ) By neglecting terms of order δx2and higher and by setting F(x+δx)=0,w e obtainaset oflinearequationsforthecorrections δxthatmoveeach functioncloser to zero simultaneously, namely J·δx=−F. (9.6.6 ) Matrix equation (9.6.6) can be solved by LUdecomposition as described in §2.3. The corrections are then added to the solution vector, xnew =xold+δx (9.6.7 ) and the process is iterated to convergence. In general it is a good idea to check the degree to which both functions and variables have converged. Once either reaches machine accuracy, the other won ’t change. Thefollowingroutine mnewtperforms ntrialiterationsstartingfromaninitial guess at the solutionvector xoflength nvariables. Iterationstops if eitherthe sum ofthemagnitudesofthefunctions Fiislessthansometolerance tolf,orthesumof theabsolutevaluesofthecorrectionsto δx iislessthansometolerance tolx.mnewt callsausersuppliedsubroutine usrfunwhichmustreturnthefunctionvalues Fand the Jacobian matrix J.I fJis difficult to compute analytically, you can try having usrfuncall the routine fdjacof§9.7 to compute the partial derivatives by finite differences. You should not make ntrialtoo big; rather inspect to see what is happening before continuing for some further iterations. SUBROUTINE mnewt(ntrial,x,n,tolx,tolf) INTEGER n,ntrial,NPREAL tolf,tolx,x(n)PARAMETER (NP=15) Up to NPvariables. C USES lubksb,ludcmp,usrfun Given an initial guess xfor a root in ndimensions, take ntrial Newton-Raphson steps to improve the root. Stop if the root converges in either summed absolute variable increments tolx or summed absolute function values tolf . INTEGER i,k,indx(NP)REAL d,errf,errx,fjac(NP,NP),fvec(NP),p(NP)do 14k=1,ntrial call usrfun(x,n,NP,fvec,fjac) User subroutine supplies function values at xinfvec and Jacobian matrix in fjac . errf=0. do11i=1,n Check function convergence. errf=errf+abs(fvec(i)) enddo 11 if(errf.le.tolf)returndo 12i=1,n Right-hand side of linear equations. p(i)=-fvec(i) enddo 12 call ludcmp(fjac,n,NP,indx,d) Solve linear equations using LU decomposition. call lubksb(fjac,n,NP,indx,p) 9.6Newton-RaphsonMethodforNonlinearSystems ofEquations 375Sample 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).errx=0. Check root convergence. do13i=1,n Update solution. errx=errx+abs(p(i)) x(i)=x(i)+p(i) enddo 13 if(errx.le.tolx)return enddo 14 returnEND Newton’sMethod versus Minimization In the next chapter, we will find that there areefficient general techniques for finding a minimum of a function of many variables. Why is that task (relatively) easy, while multidimensional root finding is often quite hard? Isn ’t minimization equivalentto findingazeroofan N-dimensionalgradientvector,notsodifferentfrom zeroingan N-dimensionalfunction? No! Thecomponentsofagradientvectorarenot independent,arbitraryfunctions. Rather,theyobeyso-calledintegrabilityconditionsthat are highly restrictive. Put crudely, you can always find a minimum by sliding downhill on a single surface. The test of “downhillness ”is thus one-dimensional. There is no analogous conceptual procedure for finding a multidimensional root, where“downhill”mustmeansimultaneouslydownhillin Nseparatefunctionspaces, thus allowing a multitude of trade-offs, as to how much progress in one dimension is worth compared with progress in another. It might occur to you to carry out multidimensional root finding by collapsing allthesedimensionsintoone: AddupthesumsofsquaresoftheindividualfunctionsF ito get a master function Fwhich (i) is positive de finite, and (ii) has a global minimum of zero exactly at all solutions of the original set of nonlinear equations. Unfortunately,asyouwillseeinthenextchapter,theef ficientalgorithmsfor finding minima come to rest on global and local minima indiscriminately. You will often find, to your great dissatisfaction, that your function Fhas a great number of local minima. InFigure9.6.1,forexample,thereislikelytobealocalminimumwherever thezerocontoursof fandgmakea closeapproachto eachother. Thepointlabeled Mis such a point, and one sees that there are no nearby roots. However, we will now see that sophisticated strategies for multidimensional rootfinding can in fact make use of the idea of minimizing a master function F,b y combining it with Newton ’s method applied to the full set of functions Fi. While such methods can still occasionally fail by coming to rest on a local minimum of F, they often succeed where a direct attack via Newton ’s method alone fails. The next section deals with these methods. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 14. [1] Ostrowski, A.M. 1966, Solutions of Equations and Systems of Equations , 2nd ed. (New York: Academic Press). Ortega, J., and Rheinboldt, W. 1970, Iterative Solution of Nonlinear Equations in Several Vari- ables(New York: Academic Press). 376 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).9.7 GloballyConvergent Methods for Nonlinear Systems of Equations We have seen that Newton ’s method for solving nonlinear equations has an unfortunate tendency to wander off into the wild blue yonder if the initial guess isnotsufficientlyclosetotheroot. A globalmethodisonethatconvergestoasolution from almost any starting point. In this section we will develop an algorithm that combinestherapidlocalconvergenceofNewton ’smethodwithagloballyconvergent strategy that will guarantee some progress towards the solution at each iteration. Thealgorithmis closelyrelatedtothequasi-Newtonmethodofminimizationwhichwe will describe in §10.7. Recall our discussion of §9.6: the Newton step for the set of equations F(x)=0 ( 9.7.1 ) is x new =xold+δx (9.7.2 ) where δx=−J−1·F (9.7.3 ) HereJistheJacobianmatrix. HowdowedecidewhethertoaccepttheNewtonstep δx? A reasonable strategy is to require that the step decrease |F|2=F·F. This is the same requirement we would impose if we were trying to minimize f=1 2F·F (9.7.4 ) (The1 2is for later convenience.) Every solution to (9.7.1) minimizes (9.7.4), but there may be local minima of (9.7.4) that are not solutions to (9.7.1). Thus, asalready mentioned, simply applying one of our minimum finding algorithms from Chapter 10 to (9.7.4) is nota good idea. To develop a better strategy, note that the Newton step (9.7.3) is a descent direction forf: ∇f·δx=(F·J)·(−J −1·F)=−F·F<0( 9.7.5 ) Thus our strategy is quite simple: We always first try the full Newton step, becauseoncewearecloseenoughtothesolutionwewillgetquadraticconvergence. However, we check at each iteration that the proposed step reduces f. If not, we backtrack alongtheNewtondirectionuntilwe haveanacceptablestep. Becausethe Newtonstepisadescentdirectionfor f,weareguaranteedto findanacceptablestep bybacktracking. We will discuss the backtrackingalgorithmin moredetail below. Notethatthismethodessentiallyminimizes fbytakingNewtonstepsdesigned tobringFtozero. Thisis notequivalenttominimizing fdirectlybytakingNewton steps designed to bring ∇fto zero. While the method can still occasionally fail by landing on a local minimum of f, this is quite rare in practice. The routine newt below will warn youif this happens. The remedyis to try a new starting point.