f9-4

9.4Newton-RaphsonMethodUsingDerivative 355Sample 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 if(p.gt.0.) q=-q Check whether in bounds. p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d Accept interpolation. d=p/q else d=xm Interpolation failed, use bisection. e=d endif else Bounds decreasing too slowly, use bisection. d=xme=d endif a=b Move last best guess to a. fa=fb if(abs(d) .gt. tol1) then Evaluate new trial root. b=b+d else b=b+sign(tol1,xm) endif fb=func(b) enddo 11 pause ’zbrent exceeding maximum iterations’zbrent=breturnEND CITED REFERENCES AND FURTHER READING: Brent,R.P.1973, AlgorithmsforMinimizationwithoutDerivatives (EnglewoodCliffs,NJ:Prentice- Hall), Chapters 3, 4. [1] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §7.2. 9.4 Newton-Raphson Method Using Derivative Perhapsthemostcelebratedofallone-dimensionalroot-findingroutinesis New- ton’smethod ,alsocalledthe Newton-Raphsonmethod . Thismethodisdistinguished from the methods of previous sections by the fact that it requires the evaluation of both the function f(x),andthe derivative f/prime(x), at arbitrary points x. The Newton-Raphson formula consists geometrically of extending the tangent line at a currentpoint xiuntil it crosses zero,thensetting thenextguess xi+1to theabscissa of that zero-crossing(see Figure 9.4.1). Algebraically, the method derives from thefamiliar Taylor series expansionof a functionin the neighborhoodof a point, f(x+δ)≈f(x)+f /prime(x)δ+f/prime/prime(x) 2δ2+.... (9.4.1 ) For small enough values of δ, and for well-behaved functions, the terms beyond linear are unimportant, hence f(x+δ)=0implies δ=−f(x) f/prime(x). (9.4.2 ) 356 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).Newton-Raphson is not restricted to one dimension. The method readily generalizes to multiple dimensions, as we shall see in §9.6 and §9.7, below. Far from a root, where the higher-order terms in the series areimportant, the Newton-Raphsonformulacangivegrosslyinaccurate,meaninglesscorrections. For instance, the initial guess for the root might be so far from the true root as to letthe search interval include a local maximum or minimum of the function. This can be death to the method (see Figure 9.4.2). If an iteration places a trial guess near such a local extremum, so that the first derivative nearly vanishes, then Newton- Raphson sends its solution off to limbo, with vanishingly small hope of recovery. Likemostpowerfultools,Newton-Raphsoncanbedestructiveusedininappropriatecircumstances. Figure 9.4.3 demonstrates another possible pathology. Why do we call Newton-Raphson powerful? The answer lies in its rate of convergence: Within a small distance /epsilon1ofxthe function and its derivative are approximately: f(x+/epsilon1)=f(x)+/epsilon1f /prime(x)+/epsilon12f/prime/prime(x) 2+···, f/prime(x+/epsilon1)=f/prime(x)+/epsilon1f/prime/prime(x)+···(9.4.3 ) By the Newton-Raphson formula, xi+1=xi−f(xi) f/prime(xi), (9.4.4 ) so that /epsilon1i+1=/epsilon1i−f(xi) f/prime(xi). (9.4.5 ) Whenatrialsolution xidiffersfromthetruerootby /epsilon1i,wecanuse(9.4.3)toexpress f(xi),f/prime(xi)in (9.4.4)in terms of /epsilon1iand derivativesat the root itself. The result is a recurrence relation for the deviations of the trial solutions /epsilon1i+1=−/epsilon12 if/prime/prime(x) 2f/prime(x). (9.4.6 ) Equation (9.4.6)says that Newton-Raphsonconverges quadratically (cf. equa- tion 9.2.3). Near a root, the number of significant digits approximately doubles with each step. This very strong convergencepropertymakes Newton-Raphsonthe methodofchoiceforanyfunctionwhosederivativecanbeevaluatedefficiently,and whose derivative is continuous and nonzero in the neighborhoodof a root. Even where Newton-Raphson is rejected for the early stages of convergence (because of its poor global convergence properties), it is very common to “polish up” a root with one or two steps of Newton-Raphson, which can multiply by two or four its number of significant figures! For an efficient realization of Newton-Raphsonthe user providesa routinethat evaluatesboth f(x)anditsfirstderivative f/prime(x)atthepoint x. TheNewton-Raphson formula can also be applied using a numerical difference to approximate the true local derivative, f/prime(x)≈f(x+dx)−f(x) dx. (9.4.7 ) 9.4Newton-RaphsonMethodUsingDerivative 357Sample 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 2 3 xf(x) Figure 9.4.1. Newton ’s method extrapolates the local derivative to find the next estimate of the root. In this example it works well and converges quadratically. f(x) x123 Figure 9.4.2. Unfortunate case where Newton ’s method encounters a local extremum and shoots off to outer space. Here bracketing bounds, as in rtsafe, would save the day. 358 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).xf(x) 21 Figure 9.4.3. Unfortunate case where Newton ’s method enters a nonconvergent cycle. This behavior is often encountered when the function fis obtained, in whole or in part, by table interpolation. With a better initial guess, the method would have succeeded. This is not, however, a recommended procedure for the following reasons: (i) You are doing two function evaluations per step, so at bestthe superlinear order of convergence will be only√ 2. (ii) If you take dxtoo small you will be wiped out by roundoff, while if you take it too large your order of convergence will be only linear, no better than using the initialevaluation f/prime(x0)for all subsequent steps. Therefore,Newton-Raphsonwithnumericalderivativesis(inonedimension)always dominated by the secant method of §9.2. (In multidimensions, where there is a paucity of available methods, Newton-Raphson with numerical derivatives must be taken more seriously. See §§9.6–9.7.) The following subroutine calls a user supplied subroutine funcd(x,fn,df) which returns the function value as fnand the derivative as df. We have included inputboundsontherootsimplytobeconsistentwithpreviousroot- findingroutines: Newton does not adjust bounds, and works only on local information at the point x. The bounds are used only to pick the midpoint as the first guess, and to reject the solution if it wanders outside of the bounds. FUNCTION rtnewt(funcd,x1,x2,xacc) INTEGER JMAXREAL rtnewt,x1,x2,xacc EXTERNAL funcd PARAMETER (JMAX=20) Set to maximum number of iterations. Using the Newton-Raphson method, find the root of a function known to lie in the interval[ x1,x2].T h e r o o t rtnewt will be refined until its accuracy is known within ±xacc .funcd is a user-supplied subroutine that returns both the function value and the first derivative of the function at the point x. INTEGER j REAL df,dx,f 9.4Newton-RaphsonMethodUsingDerivative 359Sample 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).rtnewt=.5*(x1+x2) Initial guess. do11j=1,JMAX call funcd(rtnewt,f,df) dx=f/dfrtnewt=rtnewt-dx if((x1-rtnewt)*(rtnewt-x2).lt.0.) * pause ’rtnewt jumped out of brackets’ if(abs(dx).lt.xacc) return Convergence. enddo 11 pause ’rtnewt exceeded maximum iterations’END While Newton-Raphson ’s global convergence properties are poor, it is fairly easytodesignafail-saferoutinethatutilizesacombinationofbisectionandNewton- Raphson. The hybrid algorithm takes a bisection step whenever Newton-Raphson wouldtakethesolutionoutofbounds,orwheneverNewton-Raphsonisnotreducingthe size of the brackets rapidly enough. FUNCTION rtsafe(funcd,x1,x2,xacc) INTEGER MAXITREAL rtsafe,x1,x2,xacc EXTERNAL funcd PARAMETER (MAXIT=100) Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, find the root of a function bracketedbetween x1andx2. The root, returned as the function value rtsafe , will be refined until its accuracy is known within ±xacc .funcd is a user-supplied subroutine which returns both the function value and the first derivative of the function. INTEGER j REAL df,dx,dxold,f,fh,fl,temp,xh,xl call funcd(x1,fl,df)call funcd(x2,fh,df)if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.)) * pause ’root must be bracketed in rtsafe’ if(fl.eq.0.)then rtsafe=x1 return else if(fh.eq.0.)then rtsafe=x2return else if(fl.lt.0.)then Orient the search so that f(xl)<0. xl=x1xh=x2 else xh=x1 xl=x2 endif rtsafe=.5*(x1+x2) Initialize the guess for root, dxold=abs(x2-x1) the “stepsize before last,” dx=dxold and the last step. call funcd(rtsafe,f,df) do 11j=1,MAXIT Loop over allowed iterations. if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0. Bisect if Newton out of range, * .or. abs(2.*f).gt.abs(dxold*df) ) then or not decreasing fast enough. dxold=dx dx=0.5*(xh-xl)rtsafe=xl+dxif(xl.eq.rtsafe)return Change in root is negligible. else Newton step acceptable. Take it. dxold=dxdx=f/df temp=rtsafe 360 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).rtsafe=rtsafe-dx if(temp.eq.rtsafe)return endif if(abs(dx).lt.xacc) return Convergence criterion. call funcd(rtsafe,f,df) The one new function evaluation per iteration. if(f.lt.0.) then Maintain the bracket on the root. xl=rtsafe else xh=rtsafe endif enddo 11 pause ’rtsafe exceeding maximum iterations’return END For many functions the derivative f/prime(x)often converges to machine accuracy beforethefunction f(x)itselfdoes. Whenthatisthecaseoneneednotsubsequently update f/prime(x). This shortcut is recommendedonlywhen you con fidentlyunderstand thegenericbehaviorofyourfunction,butitspeedscomputationswhenthederivative calculationislaborious. (Formallythismakestheconvergenceonlylinear,butifthe derivative isn ’t changing anyway, you can do no better.) Newton-Raphsonand Fractals An interesting sidelight to our repeated warnings about Newton-Raphson ’s unpredictable global convergence properties —its very rapid local convergence notwithstanding —is to investigate,for some particular equation,the set of starting values from which the method does, or doesn ’t converge to a root. Consider the simple equation z3−1=0 ( 9.4.8 ) whose single real root is z=1, but which also has complex roots at the other two cube roots of unity, exp(±2πi/ 3). Newton ’s method gives the iteration zj+1=zj−z3 j−1 3z2 j(9.4.9 ) Up to now, we have applied an iteration like equation (9.4.9) only for real starting values z0, but in fact all of the equations in this section also apply in the complexplane. Wecanthereforemapoutthecomplexplaneintoregionsfromwhich a starting value z0, iterated in equation (9.4.9), will, or won ’t, converge to z=1. Naively, we might expect to find a“basin of convergence ”somehow surrounding the root z=1. We surely do not expect the basin of convergence to fill the whole plane, because the plane must also contain regions that convergeto each of the two complex roots. In fact, by symmetry, the three regions must have identical shapes. Perhapstheywill bethreesymmetric 120◦wedges,withonerootcenteredineach? Now take a look at Figure 9.4.4, which shows the result of a numerical exploration. Thebasinofconvergencedoesindeedcover 1/3theareaofthecomplex plane, but its boundary is highly irregular —in fact,fractal. (A fractal, so called, has self-similar structurethat repeats on all scales of magni fication.) How does this 9.4Newton-RaphsonMethodUsingDerivative 361Sample 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).Figure 9.4.4. The complex zplane with real and imaginary components in the range (−2,2). The black region isthe setofpoints fromwhich Newton ’s methodconverges to the root z=1ofthe equation z3−1=0. Its shape is fractal. fractal emerge from something as simple as Newton ’s method, and an equation as simpleas(9.4.8)? TheanswerisalreadyimplicitinFigure9.4.2,whichshowedhow, on the real line, a local extremum causes Newton ’s method to shoot off to in finity. Suppose one is slightlyremoved from such a point. Then one might be shot off not to infinity, but—by luck—right into the basin of convergence of the desired root. But that means that in the neighborhoodof an extremum there must be a tiny,perhapsdistorted,copyofthebasinofconvergence —akindof“one-bounceaway ” copy. Similar logic shows that there can be “two-bounce ”copies,“three-bounce ” copies, and so on. A fractal thus emerges. Notice that, for equation(9.4.8),almost the whole real axis is in the domainof convergence for the root z=1. We say “almost”because of the peculiar discrete points on the negative real axis whose convergence is indeterminate (see figure). What happensif you start Newton ’s methodfrom one of these points? (Try it.) CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 2. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §8.4. Ortega, J., and Rheinboldt, W. 1970, Iterative Solution of Nonlinear Equations in Several Vari- ables(New York: Academic Press). Mandelbrot, B.B. 1983, The Fractal Geometry of Nature (San Francisco: W.H. Freeman). 362 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).Peitgen, H.-O., andSaupe, D. (eds.) 1988, TheScience of Fractal Images (NewYork: Springer- Verlag). 9.5 Roots of Polynomials Here we present a few methods for finding roots of polynomials. These will serve for most practical problemsinvolvingpolynomialsof low-to-moderatedegreeor for well-conditionedpolynomials of higher degree. Not as well appreciated as it ought to be is the fact that some polynomials are exceedingly ill-conditioned. The tiniest changes in a polynomial ’s coefficients can, in the worst case, send its roots sprawling all over the complex plane. (An infamous example due to Wilkinson is detailed by Acton [1].) Recall that a polynomial of degree nwill have nroots. The roots can be real or complex,and they might not be distinct. If the coef ficients of the polynomialare real, then complex roots will occur in pairs that are conjugate, i.e., if x1=a+bi is a root then x2=a−biwill also be a root. When the coef ficients are complex, the complex roots need not be related. Multipleroots,orcloselyspacedroots,producethemostdif ficultyfornumerical algorithms(see Figure9.5.1). Forexample, P(x)=( x−a)2has adoublereal root atx=a. However,we cannotbrackettherootbytheusualtechniqueofidentifying neighborhoods where the function changes sign, nor will slope-following methods such as Newton-Raphson work well, because both the function and its derivative vanish at a multiple root. Newton-Raphson maywork, but slowly, since large roundoff errors can occur. When a root is known in advance to be multiple, then special methods of attack are readily devised. Problems arise when (as is generally the case) we do not know in advance what pathology a root will display. Deflation ofPolynomials When seeking several or all roots of a polynomial, the total effort can be significantlyreducedbytheuseof deflation. Aseachroot risfound,thepolynomial is factored into a product involving the root and a reduced polynomial of degree one less than the original, i.e., P(x)=( x−r)Q(x). Since the roots of Qare exactly the remaining roots of P, the effort of finding additional roots decreases, becausewe workwith polynomialsoflower andlowerdegreeas we findsuccessive roots. Even more important, with de flation we can avoid the blunder of having our iterativemethodconvergetwicetothesame(nonmultiple)rootinsteadofseparately to two different roots. Deflation, which amounts to synthetic division, is a simple operation that acts onthearrayofpolynomialcoef ficients. Theconcisecodeforsyntheticdivisionbya monomial factor was given in §5.3 above. You can de flate complex roots either by convertingthatcodetocomplexdatatype,orelse —inthecaseofapolynomialwith real coefficientsbut possibly complexroots —by deflatingby a quadraticfactor, [x−(a+ib)] [x−(a−ib)] = x2−2ax+(a2+b2)( 9.5.1 )