f9-5

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 coefficients 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 coefficients are complex, the complex roots need not be related. Multipleroots,orcloselyspacedroots,producethemostdifficultyfornumerical 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 deflation we can avoid the blunder of having ouriterativemethodconvergetwicetothesame(nonmultiple)rootinsteadofseparately to two different roots. Deflation, which amounts to synthetic division, is a simple operation that acts onthearrayofpolynomialcoefficients. Theconcisecodeforsyntheticdivisionbya monomial factor was given in §5.3 above. You can deflate complex roots either by convertingthatcodetocomplexdatatype,orelse—inthecaseofapolynomialwith real coefficientsbut possibly complexroots — by deflatingby a quadraticfactor, [x−(a+ib)] [x−(a−ib)] = x 2−2ax+(a2+b2)( 9.5.1 ) 9.5RootsofPolynomials 363Sample 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)x x (b)f(x) f(x) Figure 9.5.1. (a) Linear, quadratic, and cubic behavior at the roots of polynomials. Only under high magnification (b) does it become apparent that the cubic has one, not three, roots, and that the quadratic has two roots rather than none. The routine poldivin§5.3 can be used to dividethe polynomialby this factor. Deflationmust,however,beutilizedwithcare. Becauseeachnewrootisknown with only finite accuracy, errors creep into the determination of the coef ficients of thesuccessivelyde flatedpolynomial. Consequently,therootscanbecomemoreand more inaccurate. It matters a lot whether the inaccuracy creeps in stably (plus or minusafewmultiplesofthemachineprecisionateachstage)orunstably(erosionofsuccessivesigni ficantfiguresuntiltheresultsbecomemeaningless). Whichbehavior occurs depends on just how the root is divided out. Forward de flation, where the new polynomial coef ficients are computedin the order from the highest power of x down to the constant term, was illustrated in §5.3. This turns out to be stable if the rootofsmallestabsolutevalueisdividedoutateachstage. Alternatively,onecando backwardde flation,wherenewcoef ficientsarecomputedin orderfromtheconstant term up to the coef ficient of the highest power of x. This is stable if the remaining root oflargestabsolute value is divided out at each stage. A polynomial whose coef ficients are interchanged “end-to-end, ”so that the constant becomes the highest coef ficient, etc., has its roots mapped into their reciprocals. (Proof: Divide the whole polynomial by its highest power x nand rewriteitasapolynomialin 1/x.) Thealgorithmforbackwardde flationistherefore virtuallyidenticaltothatofforwardde flation,exceptthattheoriginalcoef ficientsare taken in reverse orderand the reciprocal of the de flating root is used. Since we will useforwardde flationbelow,weleavetoyoutheexerciseofwritingaconcisecoding forbackwardde flation(asin §5.3). Formoreonthestabilityofde flation,consult [2]. To minimize the impact of increasing errors (even stable ones) when using deflation, it is advisable to treat roots of the successively de flated polynomials as onlytentativerootsoftheoriginalpolynomial. Onethen polishesthesetentativeroots by taking them as initial guesses that are to be re-solved for, using the nondeflated originalpolynomial P. Againyoumustbewarelesttwode flatedrootsareinaccurate enoughthat,underpolishing,theybothconvergetothesameunde flatedroot;inthat caseyougainaspuriousroot-multiplicityandloseadistinctroot. Thisisdetectable, 364 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).sinceyoucancompareeachpolishedrootforequalitytopreviousonesfromdistinct tentative roots. When it happens, you are advised to de flate the polynomial just once(andforthis rootonly),then againpolishthe tentativeroot,orto use Maehly ’s procedure (see equation 9.5.29 below). Belowwe saymoreabouttechniquesforpolishingrealandcomplex-conjugate tentative roots. First, let ’s get back to overall strategy. There are two schools of thought about how to proceed when faced with a polynomial of real coef ficients. One school says to go after the easiest quarry, the real,distinctroots,bythesamekindsofmethodsthatwehavediscussedinprevious sections for general functions, i.e., trial-and-error bracketing followed by a safeNewton-Raphsonas in rtsafe. Sometimes you are onlyinterested in real roots, in which case the strategy is complete. Otherwise, you then go after quadratic factors of the form (9.5.1) by any of a variety of methods. One such is Bairstow ’s method, which we will discuss below in the context of root polishing. Another is Muller ’s method, which we here brie fly discuss. Muller’sMethod Muller’smethodgeneralizesthesecantmethod,butusesquadraticinterpolation among three points instead of linear interpolation between two. Solving for thezerosof thequadraticallows the methodto findcomplexpairs ofroots. Given three previousguesses for the root x i−2,xi−1,xi, and the values ofthe polynomial P(x) at thosepoints,the nextapproximation xi+1is producedbythefollowingformulas, q≡xi−xi−1 xi−1−xi−2 A≡qP(xi)−q(1 + q)P(xi−1)+ q2P(xi−2) B≡(2q+1 ) P(xi)−(1 + q)2P(xi−1)+ q2P(xi−2) C≡(1 + q)P(xi)(9.5.2 ) followed by xi+1=xi−(xi−xi−1)/bracketleftbigg2C B±√ B2−4AC/bracketrightbigg (9.5.3 ) where the sign in the denominator is chosen to make its absolute value or modulus as large as possible. You can start the iterations with any three values of xthat you like, e.g., three equally spaced values on the real axis. Note that you must allow for the possibility of a complex denominator, and subsequent complex arithmetic, in implementing the method. Muller’s method is sometimes also used for finding complex zeros of analytic functions (not just polynomials) in the complex plane, for example in the IMSL routine ZANLY[3]. 9.5RootsofPolynomials 365Sample 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).Laguerre’s Method Thesecondschoolregardingoverallstrategyhappenstobetheonetowhichwe belong. That school advises you to use one of a verysmall numberof methods that will converge (though with greater or lesser ef ficiency) to all types of roots: real, complex, single, or multiple. Use such a method to get tentative values for all n rootsof your nth degreepolynomial. Thengobackandpolishthem as youdesire. Laguerre’smethod is byfarthemoststraightforwardofthesegeneral,complex methods. It does require complex arithmetic, even while converging to real roots; however, for polynomials with all real roots, it is guaranteed to converge to aroot from any starting point. For polynomials with some complex roots, little is theoretically proved about the method ’s convergence. Much empirical experience, however,suggeststhatnonconvergenceisextremelyunusual,and,further,canalmostalways be fixed by a simple scheme to break a nonconverginglimit cycle. (This is implemented in our routine, below.) An example of a polynomial that requires this cycle-breakingschemeis oneofhighdegree( >∼20),withallits rootsjustoutsideof the complex unit circle, approximatelyequally spaced around it. When the method convergesonasimple complexzero,it is knownthatits convergenceis thirdorder. In some instances the complex arithmetic in the Laguerre method is no disadvantage, since the polynomial itself may have complex coef ficients. Tomotivate(althoughnotrigorouslyderive)theLaguerreformulaswecannote the following relations between the polynomial and its roots and derivatives P n(x)=( x−x1)(x−x2)...(x−xn)( 9.5.4 ) ln|Pn(x)|=l n|x−x1|+l n|x−x2|+... +l n|x−xn| (9.5.5 ) dln|Pn(x)| dx=+1 x−x1+1 x−x2+... +1 x−xn=P/prime n Pn≡G (9.5.6 ) −d2ln|Pn(x)| dx2=+1 (x−x1)2+1 (x−x2)2+... +1 (x−xn)2 =/bracketleftbiggP/prime n Pn/bracketrightbigg2 −P/prime/prime n Pn≡H (9.5.7 ) Startingfromtheserelations,theLaguerreformulasmakewhatActon [1]nicelycalls “a rather drastic set of assumptions ”: The root x1that we seek is assumed to be locatedsomedistance afromourcurrentguess x, whileallother roots are assumed to be located at a distance b x−x1=a ;x−xi=bi =2 ,3,...,n (9.5.8 ) Then we can express (9.5.6), (9.5.7) as 1 a+n−1 b=G (9.5.9 ) 1 a2+n−1 b2=H (9.5.10 ) which yields as the solution for a a=n G±/radicalbig (n−1)(nH−G2)(9.5.11 ) 366 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).where the sign should be taken to yield the largest magnitude for the denominator. Since the factor inside the square root can be negative, acan be complex. (A more rigorous justi fication of equation 9.5.11 is in [4].) Themethodoperatesiteratively: Foratrialvalue x,ais calculatedbyequation (9.5.11). Then x−abecomes the next trial value. This continues until ais sufficiently small. The following routine implements the Laguerre method to find one root of a givenpolynomialofdegree m,whosecoef ficientscanbecomplex. Asusual,the first coefficient a(1)is the constant term, while a(m+1)is the coef ficient of the highest power of x. The routine implements a simpli fied version of an elegant stopping criterion due to Adams [5], which neatly balances the desire to achieve full machine accuracy, on the one hand, with the danger of iterating forever in the presence of roundoff error, on the other. SUBROUTINE laguer(a,m,x,its) INTEGER m,its,MAXIT,MR,MT REAL EPSS COMPLEX a(m+1),xPARAMETER (EPSS=2.e-7,MR=8,MT=10,MAXIT=MT*MR) Giventhedegree mandthecomplex coefficients a(1:m+1) ofthepolynomial/summationtextm+1 i=1a(i)xi−1, and given a complex value x, this routine improves xby Laguerre’s method until it con- verges, within the achievable roundoff limit, to a root of the given polynomial. The numberof iterations taken is returned as its. Parameters: EPSSis the estimated fractional roundoff error. We try to break(rare) limit cycles with MRdifferent fractional values, once every MTsteps, for MAXITtotal allowed iterations. INTEGER iter,j REAL abx,abp,abm,err,frac(MR)COMPLEX dx,x1,b,d,f,g,h,sq,gp,gm,g2SAVE frac DATA frac /.5,.25,.75,.13,.38,.62,.88,1./ Fractions used to breaka limit cycle. do 12iter=1,MAXIT Loop over iterations up to allowed maximum. its=iter b=a(m+1) err=abs(b)d=cmplx(0.,0.)f=cmplx(0.,0.) abx=abs(x) do 11j=m,1,-1 Efficient computation of the polynomial and its first two derivatives. fstores P/prime/prime/2. f=x*f+d d=x*d+b b=x*b+a(j) err=abs(b)+abx*err enddo 11 err=EPSS*err Estimate of roundoff error in evaluating polynomial. if(abs(b).le.err) then We are on the root. return else The generic case: use Laguerre’s formula. g=d/b g2=g*gh=g2-2.*f/b sq=sqrt((m-1)*(m*h-g2)) gp=g+sqgm=g-sqabp=abs(gp) abm=abs(gm) if(abp.lt.abm) gp=gmif (max(abp,abm).gt.0.) then dx=m/gp 9.5RootsofPolynomials 367Sample 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).else dx=exp(cmplx(log(1.+abx),float(iter))) endif endifx1=x-dx if(x.eq.x1)return Converged. if (mod(iter,MT).ne.0) then x=x1 else Every so often we take a fractional step, to break any limit cycle (itself a rare occurrence). x=x-dx*frac(iter/MT) endif enddo 12 pause ’too many iterations in laguer’ Veryunusual —canoccur onlyforcomplex roots. return Try a different starting guess for the root. END Here is a driverroutinethat calls laguerin succession foreach root,performs the deflation, optionally polishes the roots by the same Laguerre method —if you are not goingto polish in some other way —andfinally sorts the roots by their real parts. (We will use this routine in Chapter 13.) SUBROUTINE zroots(a,m,roots,polish) INTEGER m,MAXMREAL EPSCOMPLEX a(m+1),roots(m) LOGICAL polish PARAMETER (EPS=1.e-6,MAXM=101) A small number and maximum anticipated value of m+1. C USES laguer Giventhedegree mandthecomplex coefficients a(1:m+1) ofthepolynomial/summationtextm+1 i=1a(i)xi−1, this routine successively calls laguer and finds all mcomplex roots. The logical variable polish should be input as .true. if polishing (also by Laguerre’s method) is desired, .false. if the roots will be subsequently polished by other means. INTEGER i,j,jj,its COMPLEX ad(MAXM),x,b,cdo 11j=1,m+1 Copy of coefficients for successive deflation. ad(j)=a(j) enddo 11 do13j=m,1,-1 Loop over each root to be found. x=cmplx(0.,0.) Start at zero to favor convergence to smallest remaining root. call laguer(ad,j,x,its) Find the root. if(abs(aimag(x)).le.2.*EPS**2*abs(real(x))) x=cmplx(real(x),0.)roots(j)=xb=ad(j+1) Forward deflation. do 12jj=j,1,-1 c=ad(jj)ad(jj)=b b=x*b+c enddo 12 enddo 13 if (polish) then do14j=1,m Polish the roots using the undeflated coefficients. call laguer(a,m,roots(j),its) enddo 14 endifdo 16j=2,m Sort roots by their real parts by straight insertion. x=roots(j)do 15i=j-1,1,-1 if(real(roots(i)).le.real(x))goto 10 roots(i+1)=roots(i) enddo 15 i=0 368 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).10 roots(i+1)=x enddo 16 returnEND Eigenvalue Methods The eigenvalues of a matrix Aare the roots of the “characteristic polynomial ” P(x)=det[A−xI]. However, as we will see in Chapter 11, root- finding is not generally an ef ficient way to find eigenvalues. Turning matters around, we can use the more ef ficient eigenvalue methods that are discussed in Chapter 11 to find the roots of arbitrary polynomials. You can easily verify (see, e.g., [6]) that the characteristic polynomial of the special m×mcompanion matrix A= −a mam+1−am−1am+1··· −a2am+1−a1am+1 10 ··· 00 01 ··· 00...... 00 ··· 10 (9.5.12 ) is equivalent to the general polynomial P(x)= m+1/summationdisplay i=1aixi−1(9.5.13 ) If the coef ficients aiare real, rather than complex,then the eigenvaluesof Acan be foundusingtheroutines balancandhqrin§§11.5–11.6(seediscussionthere). This method, implemented in the routine zrhqrfollowing, is typically about a factor 2 slowerthan zroots(above). However,forsomeclassesofpolynomials,itisamore robust technique, largely because of the fairly sophisticated convergence methods embodied in hqr. If your polynomial has real coef ficients, and you are having trouble with zroots, then zrhqris a recommended alternative. SUBROUTINE zrhqr(a,m,rtr,rti) INTEGER m,MAXM REAL a(m+1),rtr(m),rti(m)PARAMETER (MAXM=50) C USES balanc,hqr Find all the roots of a polynomial with real coefficients,/summationtextm+1 i=1a(i)xi−1, given the degree mand the coefficients a(1:m+1) . The method is to construct an upper Hessenberg matrix whose eigenvalues are the desired roots, and then use the routines balanc andhqr.T h e real and imaginary parts of the roots are returned in rtr(1:m) andrti(1:m) , respectively. INTEGER j,kREAL hess(MAXM,MAXM),xr,xi if (m.gt.MAXM.or.a(m+1).eq.0.) pause ’bad args in zrhqr’ do 12k=1,m Construct the matrix. hess(1,k)=-a(m+1-k)/a(m+1)do 11j=2,m hess(j,k)=0. enddo 11 if (k.ne.m) hess(k+1,k)=1. enddo 12 9.5RootsofPolynomials 369Sample 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).call balanc(hess,m,MAXM) Find its eigenvalues. call hqr(hess,m,MAXM,rtr,rti) do14j=2,m Sort roots by their real parts by straight insertion. xr=rtr(j)xi=rti(j) do 13k=j-1,1,-1 if(rtr(k).le.xr)goto 1rtr(k+1)=rtr(k)rti(k+1)=rti(k) enddo 13 k=0 1 rtr(k+1)=xr rti(k+1)=xi enddo 14 return END OtherSure-Fire Techniques TheJenkins-Traub method has become practically a standard in black-box polynomialroot- finders,e.g.,intheIMSLlibrary [3]. Themethodistoocomplicated to discuss here, but is detailed, with referencesto the primary literature, in [4]. TheLehmer-Schur algorithm is one of a class of methods that isolate roots in the complex plane by generalizing the notion of one-dimensional bracketing. It is possible to determine ef ficiently whether there are any polynomial roots within a circle of given center and radius. From then on it is a matter of bookkeeping to hunt down all the roots by a series of decisions regarding where to place new trialcircles. Consult [1]for an introduction. Techniques for Root-Polishing Newton-Raphson works very well for real roots once the neighborhood of a root has been identi fied. The polynomial and its derivative can be ef ficiently simultaneouslyevaluatedasin §5.3. Forapolynomialofdegree n-1withcoefficients c(1)...c(n) , the following segment of code embodies one cycle of Newton- Raphson: p=c(n)*x+c(n-1) p1=c(n)do 11i=n-2,1,-1 p1=p+p1*x p=c(i)+p*x enddo 11 if (p1.eq.0.) pause ’derivative should not vanish’x=x-p/p1 Once all real roots of a polynomial have been polished, one must polish the complex roots, either directly, or by looking for quadratic factors. DirectpolishingbyNewton-Raphsonisstraightforwardforcomplexrootsifthe above code is converted to complex data types. With real polynomial coef ficients, note that your starting guess (tentative root) mustbe off the real axis, otherwise you will never get off that axis —and may get shot off to in finity by a minimum or maximum of the polynomial. 370 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).For real polynomials, the alternative means of polishing complex roots (or, for that matter,double realroots)is Bairstow’smethod,which seeks quadraticfactors. Theadvantage of going after quadratic factors is that it avoids all complex arithmetic. Bairstow ’s method seeks a quadratic factor that embodies the two roots x=a±ib, namely x2−2ax+(a2+b2)≡x2+Bx +C (9.5.14 ) In general if we divide a polynomial by a quadratic factor, there willbe a linear remainder P(x)=( x2+Bx +C)Q(x)+Rx +S. (9.5.15 ) Given BandC,RandScan be readily found, by polynomial division ( §5.3). We can consider RandStobeadjustable functionsof BandC,and theywillbezeroifthequadratic factor is a divisor of P(x). In the neighborhood of a root a first-order Taylor series expansion approximates the variation of R, Swith respect to small changes in B,C R(B+δB,C +δC)≈R(B,C )+∂R ∂BδB +∂R ∂CδC (9.5.16 ) S(B+δB,C +δC)≈S(B,C )+∂S ∂BδB +∂S ∂CδC (9.5.17 ) Toevaluate thepartialderivatives, consider the derivative of(9.5.15) withrespect to C. Since P(x)is afixed polynomial, it is independent of C, hence 0=( x2+Bx +C)∂Q ∂C+Q(x)+∂R ∂Cx+∂S ∂C(9.5.18 ) which can be rewritten as −Q(x)=( x2+Bx +C)∂Q ∂C+∂R ∂Cx+∂S ∂C(9.5.19 ) Similarly, P(x)is independent of B, so differentiating (9.5.15) with respect to Bgives −xQ (x)=( x2+Bx +C)∂Q ∂B+∂R ∂Bx+∂S ∂B(9.5.20 ) Now note that equation (9.5.19) matches equation (9.5.15) in form. Thus if we perform a secondsyntheticdivisionof P(x),i.e.,adivisionof Q(x),yieldingaremainder R1x+S1,then ∂R ∂C=−R1∂S ∂C=−S1 (9.5.21 ) To get the remaining partial derivatives, evaluate equation (9.5.20) at the two roots of the quadratic, x+andx−. Since Q(x±)=R1x±+S1 (9.5.22 ) we get ∂R ∂Bx++∂S ∂B=−x+(R1x++S1)( 9.5.23 ) ∂R ∂Bx−+∂S ∂B=−x−(R1x−+S1)( 9.5.24 ) Solve these two equations for the partial derivatives, using x++x−=−Bx +x−=C (9.5.25 ) andfind ∂R ∂B=BR 1−S1∂S ∂B=CR 1 (9.5.26 ) Bairstow’s method now consists of using Newton-Raphson intwo dimensions (which is actually the subject of the nextsection) to find a simultaneous zero of RandS. Synthetic division is used twice per cycle to evaluate R, Sand their partial derivatives with respect to B,C. Likeone-dimensional Newton-Raphson,themethodworkswellinthevicinityofaroot pair (real or complex), but it can fail miserably when started at a random point. We thereforerecommend it only in the context of polishing tentative complex roots. 9.5RootsofPolynomials 371Sample 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 qroot(p,n,b,c,eps) INTEGER n,NMAX,ITMAX REAL b,c,eps,p(n),TINY PARAMETER (NMAX=20,ITMAX=20,TINY=1.0e-6) C USES poldiv Given coefficients p(1:n)of a polynomial of degree n-1, and trial values for the coefficients of a quadratic factor x*x+b*x+c , improve the solution until the coefficients b,cchange by less than eps. The routine poldiv §5.3 is used. Parameters: At most NMAXcoefficients, ITMAXiterations. INTEGER iter REAL delb,delc,div,r,rb,rc,s,sb,sc,d(3),q(NMAX),qq(NMAX),rem(NMAX)d(3)=1.do 11iter=1,ITMAX d(2)=b d(1)=ccall poldiv(p,n,d,3,q,rem) s=rem(1) First division r,s. r=rem(2)call poldiv(q,n-1,d,3,qq,rem)sc=-rem(1) Second division partial r,swith respect to c. rc=-rem(2) sb=-c*rcrb=sc-b*rc div=1./(sb*rc-sc*rb) Solve 2x2 equation. delb=(r*sc-s*rc)*divdelc=(-r*sb+s*rb)*divb=b+delb c=c+delc if((abs(delb).le.eps*abs(b).or.abs(b).lt.TINY) * .and.(abs(delc).le.eps*abs(c) * .or.abs(c).lt.TINY)) return Coefficients converged. enddo 11 pause ’too many iterations in qroot’END We have already remarked on the annoyance of having two tentative roots collapse to one value under polishing. You are left not knowing whether yourpolishing procedure has lost a root, or whether there isactually a double root, which was split only by roundoff errors in your previous de flation. One solution is deflate-and-repolish; but de flation is what we are trying to avoid at the polishing stage. An alternative is Maehly’sprocedure . Maehly pointed out that the derivative of the reduced polynomial P j(x)≡P(x) (x−x1)···(x−xj)(9.5.27 ) can be written as P/prime j(x)=P/prime(x) (x−x1)···(x−xj)−P(x) (x−x1)···(x−xj)j/summationdisplay i=1(x−xi)−1(9.5.28 ) Hence one step of Newton-Raphson, taking a guess xkinto a new guess xk+1, can be written as xk+1 =xk−P(xk) P/prime(xk)−P(xk)/summationtextj i=1(xk−xi)−1(9.5.29 ) 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 example of so-called zero suppression as an alternative to true de flation. 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 fand gare 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 offand g(see Figure9.6.1). Unfortunately,the functions fand ghave,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