f9-7

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,weareguaranteedtofindanacceptablestep 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. 9.7GloballyConvergentMethodsforNonlinearSystems ofEquations 377Sample 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).Line Searches and Backtracking Whenwearenotcloseenough totheminimumof f,takingthefullNewtonstep p=δx need not decrease the function; we may move too far for the quadratic approximation to bevalid. All we are guaranteed is that initially fdecreases as we move in the Newton direction. So the goal is to move to a new point x newalong the direction of the Newton step p,b u t not necessarily all the way: xnew =xold+λp, 0<λ≤1( 9.7.6 ) The aim is to find λso that f(xold +λp)has decreased sufficiently. Until the early 1970s, standardpracticewastochoose λsothatxnewexactlyminimizes finthedirection p. However, we now know that it is extremely wasteful of function evaluations to do so. A better strategyisasfollows: Since pisalwaystheNewtondirectioninouralgorithms,wefirsttry λ=1,the full Newton step. This will lead to quadratic convergence when xis sufficiently close to the solution. However, if f(x new )does not meet our acceptance criteria, we backtrack along the Newtondirection,tryingasmallervalueof λ,untilwefindasuitablepoint. SincetheNewton direction is a descent direction, we are guaranteed to decrease ffor sufficiently small λ. What should the criterion for accepting a step be? It is notsufficient to require merely thatf(xnew )<f (xold). This criterion can fail to converge to a minimum of fin one of two ways. First, it is possible to construct a sequence of steps satisfying this criterion withfdecreasing too slowly relative to the step lengths. Second, one can have a sequence where the step lengths are too small relative to the initial rate of decrease of f. (For examples of such sequences, see [1], p. 117.) A simple way to fix the first problem is to require the averagerate of decrease of fto be at least some fraction αof theinitialrate of decrease ∇f·p: f(xnew )≤f(xold)+α∇f·(xnew−xold)( 9.7.7 ) Here the parameter αsatisfies 0<α< 1. We can get away with quite small values of α;α=1 0−4is a good choice. The second problem can be fixed by requiring the rate of decrease of fatxnewto be greater than some fraction βof the rate of decrease of fatxold. In practice, we will not need to impose this second constraint because our backtracking algorithm will have a built-incutoff to avoid taking steps that are too small. Here is the strategy for a practical backtracking routine: Define g(λ)≡f(x old+λp)( 9.7.8 ) so that g/prime(λ)=∇f·p (9.7.9 ) If we need to backtrack, then we model gwith the most current information we have and choose λto minimize the model. We start with g(0)andg/prime(0)available. The first step is always the Newton step, λ=1. If this step is not acceptable, we have available g(1)as well. We can therefore model g(λ)as a quadratic: g(λ)≈[g(1)−g(0)−g/prime(0)]λ2+g/prime(0)λ+g(0) ( 9.7.10 ) Taking the derivative of this quadratic, we find that it is a minimum when λ=−g/prime(0) 2[g(1)−g(0)−g/prime(0)](9.7.11 ) Since the Newton step failed, we can show that λ<∼1 2for small α. We need to guard against too small a value of λ, however. We set λmin =0.1. On second and subsequent backtracks, we model gas a cubic in λ, using the previous value g(λ1)and the second most recent value g(λ2): g(λ)=aλ3+bλ2+g/prime(0)λ+g(0) ( 9.7.12 ) 378 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).Requiring this expression to give the correct values of gatλ1andλ2gives two equations that can be solved for the coefficients aandb: /bracketleftbigga b/bracketrightbigg =1 λ1−λ2/bracketleftbigg1/λ2 1−1/λ2 2 −λ2/λ2 1λ1/λ22/bracketrightbigg ·/bracketleftbiggg(λ1)−g/prime(0)λ1−g(0) g(λ2)−g/prime(0)λ2−g(0)/bracketrightbigg (9.7.13 ) The minimum of the cubic (9.7.12) is at λ=−b+/radicalbig b2−3ag/prime(0) 3a(9.7.14 ) We enforce that λlie between λmax =0.5λ1andλmin =0.1λ1. Theroutinehastwoadditionalfeatures,aminimumsteplength alaminandamaximum step length stpmax.lnsrchwill also be used in the quasi-Newton minimization routine dfpminin the next section. SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func) INTEGER n LOGICAL checkREAL f,fold,stpmax,g(n),p(n),x(n),xold(n),func,ALF,TOLXPARAMETER (ALF=1.e-4,TOLX=1.e-7) EXTERNAL func C USES func Given an n-dimensional point xold(1:n) , the value of the function and gradient there, fold andg(1:n) , and a direction p(1:n) , finds a new point x(1:n) along the direction pfrom xold where the function func has decreased “sufficiently.” The new function value is returned in f.stpmax is an input quantity that limits the length of the steps so that you do not try to evaluate the function in regions where it is undefined or subject to overflow. pis usually the Newton direction. The output quantity check is false on a normal exit. It is true when xis too close to xold . In a minimization algorithm, this usually signals convergence and can be ignored. However, in a zero-finding algorithm the calling program should check whether the convergence is spurious. Parameters: ALF ensures sufficient decrease in function value; TOLX is the convergence criterion on ∆x. INTEGER i REAL a,alam,alam2,alamin,b,disc,f2,rhs1,rhs2,slope, * sum,temp,test,tmplam check=.false. sum=0. do11i=1,n sum=sum+p(i)*p(i) enddo 11 sum=sqrt(sum)if(sum.gt.stpmax)then Scale if attempted step is too big. do 12i=1,n p(i)=p(i)*stpmax/sum enddo 12 endif slope=0. do13i=1,n slope=slope+g(i)*p(i) enddo 13 if(slope.ge.0.) pause ’roundoff problem in lnsrch’ test=0. Compute λmin. do14i=1,n temp=abs(p(i))/max(abs(xold(i)),1.) if(temp.gt.test)test=temp enddo 14 alamin=TOLX/testalam=1. Always try full Newton step first. 1 continue Start of iteration loop. do 15i=1,n x(i)=xold(i)+alam*p(i) enddo 15 9.7GloballyConvergentMethodsforNonlinearSystems ofEquations 379Sample 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(x) if(alam.lt.alamin)then Convergence on ∆x. For zero finding, the calling program should verify the convergence.do16i=1,n x(i)=xold(i) enddo 16 check=.true.return else if(f.le.fold+ALF*alam*slope)then Sufficient function decrease. return else Backtrack. if(alam.eq.1.)then First time. tmplam=-slope/(2.*(f-fold-slope)) else Subsequent backtracks. rhs1=f-fold-alam*slope rhs2=f2-fold-alam2*slopea=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/ * (alam-alam2) if(a.eq.0.)then tmplam=-slope/(2.*b) else disc=b*b-3.*a*slopeif(disc.lt.0.)then tmplam=.5*alam else if(b.le.0.)then tmplam=(-b+sqrt(disc))/(3.*a) else tmplam=-slope/(b+sqrt(disc)) endif endif if(tmplam.gt..5*alam)tmplam=.5*alam λ≤0.5λ 1. endif endifalam2=alam f2=f alam=max(tmplam,.1*alam) λ≥0.1λ 1. goto 1 Try again. END Here now is the globally convergent Newton routine newtthat uses lnsrch. A feature ofnewtisthatyouneednotsupplytheJacobianmatrixanalytically;theroutinewillattemptto compute thenecessary partialderivatives of Fby finitedifferences in theroutine fdjac. This routineusessomeofthetechniquesdescribedin §5.7forcomputingnumericalderivatives. Of course, you can always replace fdjacwith a routine that calculates the Jacobian analytically if this is easy for you to do. SUBROUTINE newt(x,n,check) INTEGER n,nn,NP,MAXITSLOGICAL check REAL x(n),fvec,TOLF,TOLMIN,TOLX,STPMX PARAMETER (NP=40,MAXITS=200,TOLF=1.e-4,TOLMIN=1.e-6,TOLX=1.e-7, * STPMX=100.) COMMON /newtv/ fvec(NP),nn Communicates with fmin . SAVE /newtv/ C USES fdjac,fmin,lnsrch,lubksb,ludcmp Given an initial guess x(1:n) for a root in ndimensions, find the root by a globally convergent Newton’s method. The vector of functions to be zeroed, called fvec(1:n) in the routine below, is returned by a user-supplied subroutine that must be called funcv and have the declaration subroutine funcv(n,x,fvec) . The output quantity check is false on a normal return and true if the routine has converged to a local minimum of the function fmin defined below. In this case try restarting from a different initial guess. Parameters: NPis the maximum expected value of n;MAXITS is the maximum number of iterations; TOLF sets the convergence criterion on function values; TOLMIN sets the criterion for deciding whether spurious convergence to a minimum of fmin has occurred; TOLX is 380 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).the convergence criterion on δx;STPMX is the scaled maximum step length allowed in line searches. INTEGER i,its,j,indx(NP) REAL d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP), * g(NP),p(NP),xold(NP),fmin EXTERNAL fmin nn=nf=fmin(x) The vector fvec is also computed by this call. test=0. Test for initial guess being a root. Use more strin- gent test than simply TOLF . do 11i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) enddo 11 if(test.lt..01*TOLF)then check=.false. return endif sum=0. Calculate stpmax for line searches. do12i=1,n sum=sum+x(i)**2 enddo 12 stpmax=STPMX*max(sqrt(sum),float(n))do 21its=1,MAXITS Start of iteration loop. call fdjac(n,x,fvec,NP,fjac) If analytic Jacobian is available, you can replace the routine fdjac below with your own routine. do14i=1,n Compute ∇ffor the line search. sum=0. do13j=1,n sum=sum+fjac(j,i)*fvec(j) enddo 13 g(i)=sum enddo 14 do15i=1,n Storex, xold(i)=x(i) enddo 15 fold=f andf. do16i=1,n Right-hand side for linear equations. p(i)=-fvec(i) enddo 16 call ludcmp(fjac,n,NP,indx,d) Solve linear equations by LUdecomposition. call lubksb(fjac,n,NP,indx,p) call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin) lnsrch returns new xandf. It also calculates fvec at the new xwhen it calls fmin . test=0. Test for convergence on function values. do17i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) enddo 17 if(test.lt.TOLF)then check=.false. return endifif(check)then Check for gradient of fzero, i.e., spurious con- vergence. test=0. den=max(f,.5*n)do 18i=1,n temp=abs(g(i))*max(abs(x(i)),1.)/den if(temp.gt.test)test=temp enddo 18 if(test.lt.TOLMIN)then check=.true. else check=.false. endif return 9.7GloballyConvergentMethodsforNonlinearSystems ofEquations 381Sample 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 test=0. Test for convergence on δx. do19i=1,n temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.)if(temp.gt.test)test=temp enddo 19 if(test.lt.TOLX)return enddo 21 pause ’MAXITS exceeded in newt’ END SUBROUTINE fdjac(n,x,fvec,np,df) INTEGER n,np,NMAX REAL df(np,np),fvec(n),x(n),EPSPARAMETER (NMAX=40,EPS=1.e-4) C USES funcv Computes forward-difference approximation to Jacobian. On input, x(1:n) is the point at which the Jacobian is to be evaluated, fvec(1:n) is the vector of function values at the point, and npis the physical dimension of the Jacobian array df(1:n,1:n) which is output. subroutine funcv(n,x,f) is a fixed-name, user-supplied routine that returns the vector of functions at x. Parameters: NMAX is the maximum value of n;EPS is the approximate square root of the machine precision. INTEGER i,j REAL h,temp,f(NMAX)do 12j=1,n temp=x(j) h=EPS*abs(temp) if(h.eq.0.)h=EPSx(j)=temp+h Trick to reduce finite precision error. h=x(j)-temp call funcv(n,x,f)x(j)=tempdo 11i=1,n Forward difference formula. df(i,j)=(f(i)-fvec(i))/h enddo 11 enddo 12 returnEND FUNCTION fmin(x) INTEGER n,NP REAL fmin,x(*),fvecPARAMETER (NP=40) COMMON /newtv/ fvec(NP),n SAVE /newtv/ C USES funcv Returns f=1 2F·Fatx.subroutine funcv(n,x,f) is a fixed-name, user-supplied routine that returns the vector of functions at x. The common block newtv communicates the function values back to newt . INTEGER iREAL sum call funcv(n,x,fvec) sum=0.do 11i=1,n sum=sum+fvec(i)**2 enddo 11 fmin=0.5*sumreturn END Theroutine newtassumesthattypicalvaluesofallcomponentsof xandofFareoforder unity, and it can fail if this assumption is badly violated. You should rescale the variables bytheir typical values before invoking newtif this problem occurs. 382 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).MultidimensionalSecant Methods: Broyden’s Method Newton’s method as implemented above is quite powerful, but it still has several disadvantages. One drawback is that the Jacobian matrix is needed. In many problemsanalytic derivatives are unavailable. If function evaluation is expensive, then the cost offinite-difference determination of the Jacobian can be prohibitive. Justasthequasi-Newtonmethodstobediscussedin §10.7providecheapapproximations for the Hessian matrix in minimization algorithms, there are quasi-Newton methods thatprovidecheapapproximationstotheJacobianforzerofinding. Thesemethodsareoftencalledsecantmethods ,sincetheyreducetothesecantmethod( §9.2)inonedimension(see,e.g., [1]). The best of these methods still seems to be the first one introduced, Broyden’s method [2]. Let us denote the approximate Jacobian by B. Then the ith quasi-Newton step δxi is the solution of Bi·δxi=−Fi (9.7.15 ) where δxi=xi+1−xi(cf. equation 9.7.3). The quasi-Newton or secant condition is that Bi+1satisfy Bi+1·δxi=δFi (9.7.16 ) where δFi=Fi+1−Fi. Thisisthegeneralization oftheone-dimensional secantapproxima- tion to the derivative, δF/δx. However, equation (9.7.16) does not determine Bi+1uniquely in more than one dimension. Many different auxiliary conditions to pin down Bi+1have been explored, but the best-performing algorithm in practice results from Broyden’s formula. This formula is basedon the idea of getting B i+1by making the least change to Biconsistent with the secant equation (9.7.16). Broyden showed that the resulting formula is Bi+1=Bi+(δFi−Bi·δxi)⊗δxi δxi·δxi(9.7.17 ) You can easily check that Bi+1satisfies (9.7.16). Early implementations of Broyden’s method used the Sherman-Morrison formula, equation (2.7.2), to invert equation (9.7.17) analytically, B−1 i+1=B−1 i+(δxi−B−1 i·δFi)⊗δxi·B−1 i δxi·B−1 i·δFi(9.7.18 ) Then instead of solving equation (9.7.3) by e.g., LUdecomposition, one determined δxi=−B−1 i·Fi (9.7.19 ) by matrix multiplication in O(N2)operations. The disadvantage of this method is that it cannot easily be embedded in a globally convergent strategy, for which the gradient ofequation (9.7.4) requires B, notB −1, ∇(1 2F·F)/similarequalBT·F (9.7.20 ) Accordingly, we implement the update formula in the form (9.7.17). However,wecanstillpreservethe O(N2)solutionof(9.7.3)byusing QRdecomposition (§2.10)insteadof LUdecomposition. Thereasonisthatbecauseofthespecialformofequation (9.7.17), the QRdecomposition of Bican be updated into the QRdecomposition of Bi+1in O(N2)operations ( §2.10). Allweneed isaninitialapproximation B0tostarttheballrolling. Itisoften acceptable to startsimplywiththe identitymatrix, and then allow O(N)updates to produce a reasonable approximation to the Jacobian. We prefer to spend the first Nfunction evaluations on a finite-difference approximation to initialize Bvia a call to fdjac. SinceBisnottheexactJacobian,wearenotguaranteedthat δxisadescentdirectionfor f=1 2F·F(cf.equation9.7.5). Thusthelinesearchalgorithmcanfailtoreturnasuitablestep ifBwandersfarfromthetrueJacobian. Inthiscase,wereinitialize Bbyanothercallto fdjac. Like the secant method in one dimension, Broyden’s method converges superlinearly once you get close enough to the root. Embedded in a global strategy, it is almost as robust 9.7GloballyConvergentMethodsforNonlinearSystems ofEquations 383Sample 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).as Newton’s method, and often needs far fewer function evaluations to determine a zero. Note that the final value of Bisnotalways close to the true Jacobian at the root, even when the method converges. The routine broydngiven below is very similar to newtin organization. The principal differencesaretheuseof QRdecomposition insteadof LU,andtheupdating formulainstead of directly determining the Jacobian. The remarks at the end of newtabout scaling the variables apply equally to broydn. SUBROUTINE broydn(x,n,check) INTEGER n,nn,NP,MAXITS REAL x(n),fvec,EPS,TOLF,TOLMIN,TOLX,STPMXLOGICAL check PARAMETER (NP=40,MAXITS=200,EPS=1.e-7,TOLF=1.e-4,TOLMIN=1.e-6, * TOLX=EPS,STPMX=100.) COMMON /newtv/ fvec(NP),nn Communicates with fmin . SAVE /newtv/ C USES fdjac,fmin,lnsrch,qrdcmp,qrupdt,rsolv G i v e na ni n i t i a lg u e s s x(1:n) for a root in ndimensions, find the root by Broyden’s method embedded in a globally convergent strategy. The vector of functions to be zeroed, called fvec(1:n) in the routine below, is returned by a user-supplied subroutine that must be called funcv and have the declaration subroutine funcv(n,x,fvec) . The subroutine fdjac and the function fmin from newt are used. The output quantity check is false on a normal return and true if the routine has converged to a local minimum of the function fmin or if Broyden’s method can make no further progress. In this case try restarting from a different initial guess.Parameters: NPis the maximum expected value of n;MAXITS is the maximum number of iterations; EPS is close to the machine precision; TOLF sets the convergence criterion on function values; TOLMIN sets the criterion for deciding whether spurious convergence to a minimum of fmin has occurred; TOLX is the convergence criterion on δx;STPMX is the scaled maximum step length allowed in line searches. INTEGER i,its,j,k REAL den,f,fold,stpmax,sum,temp,test,c(NP),d(NP),fvcold(NP), * g(NP),p(NP),qt(NP,NP),r(NP,NP),s(NP),t(NP),w(NP),* xold(NP),fmin LOGICAL restrt,sing,skip EXTERNAL fminnn=n f=fmin(x) The vector fvec is also computed by this call. test=0. Test for initial guess being a root. Use more strin- gent test than simply TOLF . do 11i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) enddo 11 if(test.lt..01*TOLF)then check=.false.return endif sum=0. Calculate stpmax for line searches. do 12i=1,n sum=sum+x(i)**2 enddo 12 stpmax=STPMX*max(sqrt(sum),float(n))restrt=.true. Ensure initial Jacobian gets computed. do 42its=1,MAXITS Start of iteration loop. if(restrt)then call fdjac(n,x,fvec,NP,r) Initialize or reinitialize Jacobian in r. call qrdcmp(r,n,NP,c,d,sing) QRdecomposition of Jacobian. if(sing) pause ’singular Jacobian in broydn’do 14i=1,n FormQTexplicitly. do13j=1,n qt(i,j)=0. enddo 13 qt(i,i)=1. enddo 14 384 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).do18k=1,n-1 if(c(k).ne.0.)then do17j=1,n sum=0.do 15i=k,n sum=sum+r(i,k)*qt(i,j) enddo 15 sum=sum/c(k)do 16i=k,n qt(i,j)=qt(i,j)-sum*r(i,k) enddo 16 enddo 17 endif enddo 18 do21i=1,n FormRexplicitly. r(i,i)=d(i) do19j=1,i-1 r(i,j)=0. enddo 19 enddo 21 else Carry out Broyden update. do22i=1,n s=δx. s(i)=x(i)-xold(i) enddo 22 do24i=1,n t=R·s. sum=0.do 23j=i,n sum=sum+r(i,j)*s(j) enddo 23 t(i)=sum enddo 24 skip=.true. do26i=1,n w=δF−B·s. sum=0. do25j=1,n sum=sum+qt(j,i)*t(j) enddo 25 w(i)=fvec(i)-fvcold(i)-sum if(abs(w(i)).ge.EPS*(abs(fvec(i))+abs(fvcold(i))))then Don’t update with noisy components of w. skip=.false. else w(i)=0. endif enddo 26 if(.not.skip)then do28i=1,n t=QT·w. sum=0. do27j=1,n sum=sum+qt(i,j)*w(j) enddo 27 t(i)=sum enddo 28 den=0. do29i=1,n den=den+s(i)**2 enddo 29 do31i=1,n Stores/(s·s)ins. s(i)=s(i)/den enddo 31 call qrupdt(r,qt,n,NP,t,s) UpdateRandQT. do32i=1,n if(r(i,i).eq.0.) pause ’r singular in broydn’ d(i)=r(i,i) Diagonal of Rstored in d. 9.7GloballyConvergentMethodsforNonlinearSystems ofEquations 385Sample 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).enddo 32 endif endif do34i=1,n Right-hand side for linear equations is −QT·F. sum=0. do33j=1,n sum=sum+qt(i,j)*fvec(j) enddo 33 p(i)=-sum enddo 34 do36i=n,1,-1 Compute ∇f≈(Q·R)T·Ffor the line search. sum=0.do 35j=1,i sum=sum-r(j,i)*p(j) enddo 35 g(i)=sum enddo 36 do37i=1,n StorexandF. xold(i)=x(i)fvcold(i)=fvec(i) enddo 37 fold=f Store f. call rsolv(r,n,NP,d,p) Solve linear equations. call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin) lnsrch returns new xandf. It also calculates fvec at the new xwhen it calls fmin . test=0. Test for convergence on function values. do38i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) enddo 38 if(test.lt.TOLF)then check=.false. return endifif(check)then True if line search failed to find a new x. if(restrt)then Failure; already tried reinitializing the Jacobian. return else Check for gradient of fzero, i.e., spurious con- vergence. test=0. den=max(f,.5*n) do 39i=1,n temp=abs(g(i))*max(abs(x(i)),1.)/den if(temp.gt.test)test=temp enddo 39 if(test.lt.TOLMIN)then return else Try reinitializing the Jacobian. restrt=.true. endif endif else Successful step; will use Broyden update for next step. restrt=.false. test=0. Test for convergence on δx. do41i=1,n temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.)if(temp.gt.test)test=temp enddo 41 if(test.lt.TOLX)return endif enddo 42 pause ’MAXITS exceeded in broydn’END 386 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).More Advanced Implementations One of the principal ways that the methods described so far can fail is if J(in Newton’s method) or Bin (Broyden’s method) becomes singular or nearly singular, so that δxcannot be determined. If you are lucky, this situation will not occur very often in practice. Methodsdeveloped so far to deal with this problem involve monitoring the condition number of Jand perturbing Jif singularity or near singularity is detected. This is most easily implemented if the QRdecomposition is used instead of LUin Newton’s method (see [1]for details). Our personal experience is that, while such an algorithm can solve problems where Jis exactly singular and the standard Newton’s method fails, it is occasionally less robust onother problems where LUdecomposition succeeds. Clearlyimplementation details involving roundoff, underflow, etc., are important here and the last word is yet to be written. Our global strategies both for minimization and zero finding have been based on line searches. Other global algorithms, such as the hook step anddogleg step methods, are based instead on the model-trust region approach, which is related to the Levenberg-Marquardt algorithm for nonlinear least-squares ( §15.5). While somewhat more complicated than line searches, these methods have a reputation for robustness even when starting far from thedesired zero or minimum [1]. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). [1] Broyden, C.G. 1965, Mathematics of Computation , vol. 19, pp. 577–593. [2]