9.5 Roots of Polynomials 369 9.5 Roots of Polynomials Here we present a few methods for finding roots of polynomials.These will serve for most practical problems involving polynomials of low-to-moderate degree or for well-conditioned polynomials 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 n will have n roots.The roots can be real 81 or complex,and they might not be distinct.If the coefficients of the polynomial are real,then complex roots will occur in pairs that are conjugate,i.e.,if=a+bi 虽2 is a root then z2 =a-bi will also be a root.When the coefficients are complex, the complex roots need not be related. Multiple roots,or closely spaced roots,produce the most difficulty for numerical algorithms (see Figure 9.5.1).For example,P(x)=(x-a)2 has a double real root at z=a.However,we cannot bracket the root by the usual technique of identifying neighborhoods where the function changes sign,nor will slope-following methods 9 such as Newton-Raphson work well,because both the function and its derivative vanish at a multiple root.Newton-Raphson may work,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. 超% 9 Deflation of Polynomials When seeking several or all roots of a polynomial,the total effort can be 6 significantly reduced by the use of deftation.As each root r is found,the polynomial 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(z).Since the roots of are exactly the remaining roots of P,the effort of finding additional roots decreases, because we work with polynomials of lower and lower degree as we find successive 10621 roots.Even more important,with deflation we can avoid the blunder of having our iterative method converge twice to the same(nonmultiple)root instead of separately "hg2 Numerical Recipes 43106 to two different roots. Deflation,which amounts to synthetic division,is a simple operation that acts (outside on the array of polynomial coefficients.The concise code for synthetic division by a monomial factor was given in 85.3 above.You can deflate complex roots either by North converting that code to complex data type,or else-in the case of a polynomial with real coefficients but possibly complex roots-by deflating by a quadratic factor, [x-(a+b)][x-(a-b)]=x2-2ax+(a2+b2) (9.5.1) The routine poldiv in 85.3 can be used to divide the polynomial by this factor. Deflation must,however,be utilized with care.Because each new root is known with only finite accuracy,errors creep into the determination of the coefficients of the successively deflated polynomial.Consequently,the roots can become more and more inaccurate.It matters a lot whether the inaccuracy creeps in stably (plus or
9.5 Roots of Polynomials 369 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 9.5 Roots of Polynomials Here we present a few methods for finding roots of polynomials. These will serve for most practical problems involving polynomials of low-to-moderate degree or for well-conditioned polynomials 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 n will have n roots. The roots can be real or complex, and they might not be distinct. If the coefficients of the polynomial are real, then complex roots will occur in pairs that are conjugate, i.e., if x1 = a + bi is a root then x2 = a − bi will also be a root. When the coefficients are complex, the complex roots need not be related. Multiple roots, or closely spaced roots, produce the most difficulty for numerical algorithms (see Figure 9.5.1). For example, P(x)=(x − a) 2 has a double real root at x = a. However, we cannot bracket the root by the usual technique of identifying 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 may work, 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 of Polynomials When seeking several or all roots of a polynomial, the total effort can be significantly reduced by the use of deflation. As each root r is found, the polynomial 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 Q are exactly the remaining roots of P, the effort of finding additional roots decreases, because we work with polynomials of lower and lower degree as we find successive roots. Even more important, with deflation we can avoid the blunder of having our iterative method converge twice to the same (nonmultiple) root instead of separately to two different roots. Deflation, which amounts to synthetic division, is a simple operation that acts on the array of polynomial coefficients. The concise code for synthetic division by a monomial factor was given in §5.3 above. You can deflate complex roots either by converting that code to complex data type, or else — in the case of a polynomial with real coefficients but possibly complex roots — by deflating by a quadratic factor, [x − (a + ib)] [x − (a − ib)] = x2 − 2ax + (a2 + b2) (9.5.1) The routine poldiv in §5.3 can be used to divide the polynomial by this factor. Deflation must, however, be utilized with care. Because each new root is known with only finite accuracy, errors creep into the determination of the coefficients of the successively deflated polynomial. Consequently, the roots can become more and more inaccurate. It matters a lot whether the inaccuracy creeps in stably (plus or
370 Chapter 9.Root Finding and Nonlinear Sets of Equations fx) 八x) e小N 83 from NUMERICAL RECIPESI 188891992 11800 (a) (b) 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. 9 minus a few multiples of the machine precision at each stage)or unstably (erosion of America computer, successive significant figures until the results become meaningless).Which behavior 9。 occurs depends on just how the root is divided out.Forward deflation,where the new polynomial coefficients are computed in the order from the highest power of x 9 down to the constant term,was illustrated in $5.3.This turns out to be stable if the root of smallest absolute value is divided out at each stage.Alternatively.one can do backward deflation,where new coefficients are computed in order from the constant term up to the coefficient of the highest power of z.This is stable if the remaining 可 root of largest absolute value is divided out at each stage. A polynomial whose coefficients are interchanged "end-to-end,"so that the constant becomes the highest coefficient,etc.,has its roots mapped into their reciprocals.(Proof:Divide the whole polynomial by its highest power x and rewrite it as a polynomial in 1/z.)The algorithm for backward deflation is therefore Fuurggoglrion Numerical Recipes 10621 virtually identical to that of forward deflation,except that the original coefficients are taken in reverse order and the reciprocal of the deflating root is used.Since we will 43108 use forward deflation below,we leave to you the exercise of writing a concise coding for backward deflation (as in 85.3).For more on the stability of deflation.consult [2]. To minimize the impact of increasing errors(even stable ones)when using (outside deflation.it is advisable to treat roots of the successively deflated polynomials as North Software. only tentative roots ofthe original polynomial.One then polishes these tentative roots by taking them as initial guesses that are to be re-solved for,using the nondeflated original polynomial P.Again you must beware lest two deflated roots are inaccurate enough that,under polishing,they both converge to the same undeflated root;in that case you gain a spurious root-multiplicity and lose a distinct root.This is detectable. since you can compare each polished root for equality to previous ones from distinct tentative roots.When it happens,you are advised to deflate the polynomial just once (and for this root only),then again polish the tentative root,or to use Maehly's procedure (see equation 9.5.29 below). Below we say more about techniques for polishing real and complex-conjugate
370 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). (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. minus a few multiples of the machine precision at each stage) or unstably (erosion of successive significant figures until the results become meaningless). Which behavior occurs depends on just how the root is divided out. Forward deflation, where the new polynomial coefficients are computed in 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 root of smallest absolute value is divided out at each stage. Alternatively, one can do backward deflation, where new coefficients are computed in order from the constant term up to the coefficient of the highest power of x. This is stable if the remaining root of largest absolute value is divided out at each stage. A polynomial whose coefficients are interchanged “end-to-end,” so that the constant becomes the highest coefficient, etc., has its roots mapped into their reciprocals. (Proof: Divide the whole polynomial by its highest power x n and rewrite it as a polynomial in 1/x.) The algorithm for backward deflation is therefore virtually identical to that of forward deflation, except that the original coefficients are taken in reverse order and the reciprocal of the deflating root is used. Since we will use forward deflation below, we leave to you the exercise of writing a concise coding for backward deflation (as in §5.3). For more on the stability of deflation, consult [2]. To minimize the impact of increasing errors (even stable ones) when using deflation, it is advisable to treat roots of the successively deflated polynomials as only tentative roots of the original polynomial. One then polishesthese tentative roots by taking them as initial guesses that are to be re-solved for, using the nondeflated original polynomial P. Again you must beware lest two deflated roots are inaccurate enough that, under polishing, they both converge to the same undeflated root; in that case you gain a spurious root-multiplicity and lose a distinct root. This is detectable, since you can compare each polished root for equality to previous ones from distinct tentative roots. When it happens, you are advised to deflate the polynomial just once (and for this root only), then again polish the tentative root, or to use Maehly’s procedure (see equation 9.5.29 below). Below we say more about techniques for polishing real and complex-conjugate
9.5 Roots of Polynomials 371 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 coefficients.One school says to go after the easiest quarry,the real,distinct roots,by the same kinds of methods that we have discussed in previous sections for general functions,i.e.,trial-and-error bracketing followed by a safe Newton-Raphson as in rtsafe.Sometimes you are only interested 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 briefly discuss. Muller's Method Muller's method generalizes the secant method,but uses quadratic interpolation 茶 among three points instead of linear interpolation between two.Solving for the zeros of the quadratic allows the method to find complex pairs of roots.Given three previous guesses for the rooti-2,-1,i,and the values of the polynomial P() RECIPES I at those points,the next approximation is produced by the following formulas, q三 Ti-Ti-1 University Press. 工i-1-xi-2 A=qP(xi)-q(1+q)P(i-1)+q2P(zi_2) (9.5.2) B=(2g+1)P(x)-(1+q)2P(x-1)+q2P(x-2) C≡(1+q)P(x) SCIENTIFIC followed by 2C xi+1=工2-(x-t-1) (9.5.3) B士VB2-4AC 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 x that you Numerical Recipes 10.621 43106 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. (outside 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]. Laguerre's Method The second school regarding overall strategy happens to be the one to which we belong.That school advises you to use one of a very small number of methods that will converge(though with greater or lesser efficiency)to all types of roots:real, complex,single,or multiple.Use such a method to get tentative values for all n roots of your nth degree polynomial.Then go back and polish them as you desire
9.5 Roots of Polynomials 371 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 coefficients. One school says to go after the easiest quarry, the real, distinct roots, by the same kinds of methods that we have discussed in previous sections for general functions, i.e., trial-and-error bracketing followed by a safe Newton-Raphson as in rtsafe. Sometimes you are only interested 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 briefly discuss. Muller’s Method Muller’s method generalizes the secant method, but uses quadratic interpolation among three points instead of linear interpolation between two. Solving for the zeros of the quadratic allows the method to find complex pairs of roots. Given three previous guesses for the root xi−2, xi−1, xi, and the values of the polynomial P(x) at those points, the next approximation xi+1 is produced by the following formulas, 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) 2C B ± √ B2 − 4AC (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 x that 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]. Laguerre’s Method The second school regarding overall strategy happens to be the one to which we belong. That school advises you to use one of a very small number of methods that will converge (though with greater or lesser efficiency) to all types of roots: real, complex, single, or multiple. Use such a method to get tentative values for all n roots of your nth degree polynomial. Then go back and polish them as you desire.
372 Chapter 9.Root Finding and Nonlinear Sets of Equations Laguerre's method is by far the most straightforward of these general,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 a root from any starting point.For polynomials with some complex roots,little is theoretically proved about the method's convergence.Much empirical experience, however,suggests that nonconvergence is extremely unusual,and,further,can almost always be fixed by a simple scheme to break a nonconverging limit cycle.(This is implemented in our routine,below.)An example of a polynomial that requires this cycle-breaking scheme is one of high degree(20),with all its roots just outside of the complex unit circle,approximately equally spaced around it.When the method 81 converges on a simple complex zero,it is known that its convergence is third order. In some instances the complex arithmetic in the Laguerre method is no disadvantage,since the polynomial itself may have complex coefficients. To motivate(although not rigorously derive)the Laguerre formulas we can note 分 the following relations between the polynomial and its roots and derivatives Pn(x)=(x-x1)(x-x2)..(x-xn) (9.5.4) In Pn(x)=Inx-x1+Inz-z2+...+Inx-2n (9.5.5) 9 dx T-71 T-72 T-In P A2,9 d2In P() 1 1 1 dx2 三十1 -1产+(-2P++红-m $n∽ 9 (9.5.7) P Starting from these relations,the Laguerre formulas make what Acton [1]nicely calls "a rather drastic set of assumptions":The root xI that we seek is assumed to be located some distance a from our current guess x,while all other roots are assumed to be located at a distance b 9 x-x1=a x-zi=b i=2,3,...,n (9.5.8) Numerica 10.621 Then we can express (9.5.6),(9.5.7)as 43126 1,n-1 一十 -=G a b (9.5.9) 1n-1 a2+ 2=H (9.5.10) which yields as the solution for a Q= (9.5.11) G±√(n-1)(nH-G) where the sign should be taken to yield the largest magnitude for the denominator. Since the factor inside the square root can be negative,a can be complex.(A more rigorous justification of equation 9.5.11 is in [41.)
372 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Laguerre’s method is by far the most straightforward of these general, 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 a root from any starting point. For polynomials with some complex roots, little is theoretically proved about the method’s convergence. Much empirical experience, however, suggests that nonconvergenceis extremely unusual, and, further, can almost always be fixed by a simple scheme to break a nonconverging limit cycle. (This is implemented in our routine, below.) An example of a polynomial that requires this cycle-breaking scheme is one of high degree (>∼ 20), with all its roots just outside of the complex unit circle, approximately equally spaced around it. When the method converges on a simple complex zero, it is known that its convergence is third order. In some instances the complex arithmetic in the Laguerre method is no disadvantage, since the polynomial itself may have complex coefficients. To motivate (although not rigorously derive) the Laguerre formulas we can note the following relations between the polynomial and its roots and derivatives Pn(x)=(x − x1)(x − x2)...(x − xn) (9.5.4) ln |Pn(x)| = ln |x − x1| + ln |x − x2| + ... + ln |x − xn| (9.5.5) d ln |Pn(x)| dx = + 1 x − x1 + 1 x − x2 + ... + 1 x − xn = P n Pn ≡ G (9.5.6) −d2 ln |Pn(x)| dx2 = + 1 (x − x1)2 + 1 (x − x2)2 + ... + 1 (x − xn)2 = P n Pn 2 − P n Pn ≡ H (9.5.7) Starting from these relations, the Laguerre formulas make what Acton [1] nicely calls “a rather drastic set of assumptions”: The root x1 that we seek is assumed to be located some distance a from our current guess x, while all other roots are assumed to be located at a distance b x − x1 = a ; x − xi = b i = 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 ± (n − 1)(nH − G2) (9.5.11) where the sign should be taken to yield the largest magnitude for the denominator. Since the factor inside the square root can be negative, a can be complex. (A more rigorous justification of equation 9.5.11 is in [4].)
9.5 Roots of Polynomials 373 The method operates iteratively:For a trial value z.a is calculated by equation (9.5.11).Then x-a becomes the next trial value.This continues until a is sufficiently small. The following routine implements the Laguerre method to find one root of a given polynomial of degree m,whose coefficients can be complex.As usual,the first coefficient a[o]is the constant term,while a [m]is the coefficient of the highest power of x.The routine implements a simplified 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. #include 19881992 #include "complex.h" #include "nrutil.h" 1.800 #define EPSS 1.0e-7 #define MR 8 #define MT 10 from NUMERICAL RECIPESI #define MAXIT (MT*MR) Here EPSS is the estimated fractional roundoff error.We try to break (rare)limit cycles with MR different fractional values,once every MT steps,for MAXIT total allowed iterations. server (Nort void laguer(fcomplex a[],int m,fcomplex *x,int *its) THE Given the degreemand the m1complex coefficients a[.]of the polynomial America computer and given a complex value x,this routine improves x by Laguerre's method until it converges, ART within the achievable roundoff limit,to a root of the given polynomial.The number of iterations taken is returned as its. 9 Programs int iter,j; float abx,abp,abm,err; fcomplex dx,x1,b,d,f,g,h,sq,gp,gm,g2; stat1cf1 oat frac[R+1]={0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0J: Fractions used to break a limit cycle. for (iter=1;iter=0;j--)( Efficient computation of the polynomial and f=Cadd(Cmul(*x,f),d); its first two derivatives.f stores P/2. Numerical Recipes -43108 d=Cadd(Cmul(*x,d),b); b=Cadd(Cmul(*x,b),a[j]) err-Cabs(b)+abxerr; (outside err米=EPSS North Software. Estimate of roundoff error in evaluating polynomial. if (Cabs(b)0.0 Cdiv(Complex((float)m,0.0),gp) RCmul(1+abx,Complex(cos((float)iter),sin((float)iter))))); x1=Csub(*x,dx);
9.5 Roots of Polynomials 373 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The method operates iteratively: For a trial value x, a is calculated by equation (9.5.11). Then x − a becomes the next trial value. This continues until a is sufficiently small. The following routine implements the Laguerre method to find one root of a given polynomial of degree m, whose coefficients can be complex. As usual, the first coefficient a[0] is the constant term, while a[m] is the coefficient of the highest power of x. The routine implements a simplified 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. #include #include "complex.h" #include "nrutil.h" #define EPSS 1.0e-7 #define MR 8 #define MT 10 #define MAXIT (MT*MR) Here EPSS is the estimated fractional roundoff error. We try to break (rare) limit cycles with MR different fractional values, once every MT steps, for MAXIT total allowed iterations. void laguer(fcomplex a[], int m, fcomplex *x, int *its) Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial m i=0 a[i]xi, and given a complex value x, this routine improves x by Laguerre’s method until it converges, within the achievable roundoff limit, to a root of the given polynomial. The number of iterations taken is returned as its. { int iter,j; float abx,abp,abm,err; fcomplex dx,x1,b,d,f,g,h,sq,gp,gm,g2; static float frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0}; Fractions used to break a limit cycle. for (iter=1;iter=0;j--) { Efficient computation of the polynomial and its first two derivatives. f stores P f=Cadd(Cmul(*x,f),d); /2. d=Cadd(Cmul(*x,d),b); b=Cadd(Cmul(*x,b),a[j]); err=Cabs(b)+abx*err; } err *= EPSS; Estimate of roundoff error in evaluating polynomial. if (Cabs(b) 0.0 ? Cdiv(Complex((float) m,0.0),gp) : RCmul(1+abx,Complex(cos((float)iter),sin((float)iter))))); x1=Csub(*x,dx);
374 Chapter 9.Root Finding and Nonlinear Sets of Equations if (x->r =x1.r &x->i =x1.i)return; Converged. if (iter MT)*x-x1; else *x-Csub(*x,RCmul(frac[iter/MT],dx)); Every so often we take a fractional step,to break any limit cycle (itself a rare occur- rence). y nrerror("too many iterations in laguer"); Very unusual-can occur only for complex roots.Try a different starting guess for the root. return; Here is a driver routine that calls laguer in succession for each root,performs the deflation,optionally polishes the roots by the same Laguerre method-if you 18881892 are not going to polish in some other way-and finally sorts the roots by their real parts.(We will use this routine in Chapter 13.) from NUMERICAL RECIPESI #include #include "complex.h" #define EPS 2.0e-6 #define MAXM 100 A small number,and maximum anticipated value of m America computer, void zroots(fcomplex a[],int m,fcomplex roots[],int polish) ART Given the degree m and the m+1 complex coefficients a[0..m]of the polynomial a(i)xi this routine successively calls laguer and findsallcomplex roots in roots[The boolean variable polish should be input as true (1)if polishing (also by Laguerre's method) Programs is desired,false (0)if the roots will be subsequently polished by other means. void laguer(fcomplex a[],int m,fcomplex *x,int *its); int i,its,j,jji fcomplex x,b,c,ad[MAXM]; to dir for (j=0;j=1;j=-){ Loop over each root to be found. 1920 SCIENTIFIC COMPUTING(ISBN x=Complex(0.0,0.0) Start at zero to favor convergence to small- laguer(ad,j,&x,&its); est remaining root,and find the root. 1f(fabs(x.1)=0;jj--){ c=ad[ji]; Numerical Recipes 43108 ad[jj]=b; b=Cadd(Cmul(x,b),c) (outside if (polish) Software. for (j=1;j=1;1-){ if (roots[i].r <x.r)break; roots[i+1]=roots[i]; roots[i+1]=x; 2
374 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). if (x->r == x1.r && x->i == x1.i) return; Converged. if (iter % MT) *x=x1; else *x=Csub(*x,RCmul(frac[iter/MT],dx)); Every so often we take a fractional step, to break any limit cycle (itself a rare occurrence). } nrerror("too many iterations in laguer"); Very unusual — can occur only for complex roots. Try a different starting guess for the root. return; } Here is a driver routine that calls laguer in succession for each root, performs the deflation, optionally polishes the roots by the same Laguerre method — if you are not going to polish in some other way — and finally sorts the roots by their real parts. (We will use this routine in Chapter 13.) #include #include "complex.h" #define EPS 2.0e-6 #define MAXM 100 A small number, and maximum anticipated value of m. void zroots(fcomplex a[], int m, fcomplex roots[], int polish) Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial m i=0 a(i)xi, this routine successively calls laguer and finds all m complex roots in roots[1..m]. The boolean variable polish should be input as true (1) if polishing (also by Laguerre’s method) is desired, false (0) if the roots will be subsequently polished by other means. { void laguer(fcomplex a[], int m, fcomplex *x, int *its); int i,its,j,jj; fcomplex x,b,c,ad[MAXM]; for (j=0;j=1;j--) { Loop over each root to be found. x=Complex(0.0,0.0); Start at zero to favor convergence to smalllaguer(ad,j,&x,&its); est remaining root, and find the root. if (fabs(x.i) =0;jj--) { c=ad[jj]; ad[jj]=b; b=Cadd(Cmul(x,b),c); } } if (polish) for (j=1;j=1;i--) { if (roots[i].r <= x.r) break; roots[i+1]=roots[i]; } roots[i+1]=x; } }
9.5 Roots of Polynomials 375 Eigenvalue Methods The eigenvalues of a matrix A are the roots of the "characteristic polynomial" P(z)=det[A-xI].However,as we will see in Chapter 11,root-finding is not generally an efficient way to find eigenvalues.Turning matters around,we can use the more efficient eigenvalue methods that are discussed in Chapter 11 to find the roots of arbitrary polynomials.You can easily verify (see,e.g.,[61)that the characteristic polynomial of the special m x m companion matrix am-1 am-2 a1 am am am am 1 0 0 0 0 0 0 (9.5.12) 0 0 0 is equivalent to the general polynomial (9.5.13) America computer, Press. If the coefficients a;are real,rather than complex,then the eigenvalues of A can be ART found using the routines balanc and hgr in 8811.5-11.6(see discussion there).This Program method,implemented in the routine zrhgr following,is typically about a factor 2 slower than zroots (above).However,for some classes of polynomials,it is a more robust technique,largely because of the fairly sophisticated convergence methods embodied in hgr.If your polynomial has real coefficients,and you are having 6 d爱 trouble with zroots,then zrhgr is a recommended alternative. #include "nrutil.h" OF SCIENTIFIC COMPUTING(ISBN #define MAXM 50 1988-19920 void zrhgr(float a[],int m,float rtr[],float rti[]) Find all the roots of a polynomial with real coefficients.a(i)zi,given the degree m and the coefficients a [0..m].The method is to construct an upper Hessenberg matrix whose eigenvalues are the desired roots,and then use the routines balanc and hgr.The real and imaginary parts of the roots are returned in rtr[1..m]and rti[1..m],respectively. Numerical Recipes 10-521 43108 void balanc(float **a,int n); void hgr(float *a,int n,float wr[],float wi[]); (outside int j,k; Software. float **hess,xr,xi; hess=matrix(1,MAXM,1,MAXM) Amer ying of if (m MAXM II a[m]=0.0)nrerror("bad args in zrhgr"); for(k=1;k<=m;k++){ Construct the matrix. hess[1][k]-a[m-k]/a[m]; for (j=2;j<=m;j++)hess[j][k]=0.0; if (k !=m)hess[k+1][k]=1.0; balanc(hess,m); Find its eigenvalues. hgr(hess,m,rtr,rti); for(j=2;j<=n;j++){ Sort roots by their real parts by straight insertion. xr=rtr[j]; xi=rti[j];
9.5 Roots of Polynomials 375 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Eigenvalue Methods The eigenvalues of a matrix A are 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 efficient way to find eigenvalues. Turning matters around, we can use the more efficient 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 × m companion matrix A = −am−1 am −am−2 am ··· − a1 am − a0 am 1 0 ··· 0 0 0 1 ··· 0 0 . . . . . . 0 0 ··· 1 0 (9.5.12) is equivalent to the general polynomial P(x) = m i=0 aixi (9.5.13) If the coefficients ai are real, rather than complex, then the eigenvalues of A can be found using the routines balanc and hqr in §§11.5–11.6 (see discussion there). This method, implemented in the routine zrhqr following, is typically about a factor 2 slower than zroots (above). However, for some classes of polynomials, it is a more robust technique, largely because of the fairly sophisticated convergence methods embodied in hqr. If your polynomial has real coefficients, and you are having trouble with zroots, then zrhqr is a recommended alternative. #include "nrutil.h" #define MAXM 50 void zrhqr(float a[], int m, float rtr[], float rti[]) Find all the roots of a polynomial with real coefficients, m i=0 a(i)xi, given the degree m and the coefficients a[0..m]. The method is to construct an upper Hessenberg matrix whose eigenvalues are the desired roots, and then use the routines balanc and hqr. The real and imaginary parts of the roots are returned in rtr[1..m] and rti[1..m], respectively. { void balanc(float **a, int n); void hqr(float **a, int n, float wr[], float wi[]); int j,k; float **hess,xr,xi; hess=matrix(1,MAXM,1,MAXM); if (m > MAXM || a[m] == 0.0) nrerror("bad args in zrhqr"); for (k=1;k<=m;k++) { Construct the matrix. hess[1][k] = -a[m-k]/a[m]; for (j=2;j<=m;j++) hess[j][k]=0.0; if (k != m) hess[k+1][k]=1.0; } balanc(hess,m); Find its eigenvalues. hqr(hess,m,rtr,rti); for (j=2;j<=m;j++) { Sort roots by their real parts by straight insertion. xr=rtr[j]; xi=rti[j];
376 Chapter 9.Root Finding and Nonlinear Sets of Equations for(k=j-1;k>=1;k--)[ if (rtr[k]=0;1--)[ P1=p+p1*z; 之州 6 COMPUTING (ISBN 1988-19920 Numerical Recipes 10-621 43108 p=c[i]+p*x; if (p1 ==0.0)nrerror("derivative should not vanish"); x-=p/p1; (outside Software. Once all real roots of a polynomial have been polished,one must polish the complex roots,either directly,or by looking for quadratic factors Direct polishing by Newton-Raphson is straightforward for complex roots if the above code is converted to complex data types.With real polynomial coefficients, note that your starting guess (tentative root)must be off the real axis,otherwise you will never get off that axis-and may get shot off to infinity by a minimum or maximum of the polynomial. For real polynomials,the alternative means of polishing complex roots (or,for that matter,double real roots)is Bairstow's method,which seeks quadratic factors.The advantage
376 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (k=j-1;k>=1;k--) { if (rtr[k] =0;i--) { p1=p+p1*x; p=c[i]+p*x; } if (p1 == 0.0) nrerror("derivative should not vanish"); 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. Direct polishing by Newton-Raphson is straightforward for complex roots if the above code is converted to complex data types. With real polynomial coefficients, note that your starting guess (tentative root) must be off the real axis, otherwise you will never get off that axis — and may get shot off to infinity by a minimum or maximum of the polynomial. For real polynomials, the alternative means of polishing complex roots (or, for that matter, double real roots) is Bairstow’s method, which seeks quadratic factors. The advantage
9.5 Roots of Polynomials 377 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+2)=22+Bx+C (9.5.14) In general if we divide a polynomial by a quadratic factor,there will be a linear remainder P(x)=(x2+Bx+C)Q(x)+Rx+S. (9.5.15) Given B and C,R and S can be readily found,by polynomial division($5.3).We can consider R and S to be adjustable functions of B and C.and they will be zero if the quadratic factor is a divisor of P(). In the neighborhood of a root a first-order Taylor series expansion approximates the variation of R,S with respect to small changes in B,C R(B+6B.C+6C)R(B.C)+R5B+ aR8C (9.5.16) 8B ac S(B+6B.C+6C)S(B.C)+ 2iB+90 aB (9.5.17) ICAL To evaluate the partial derivatives,consider the derivative of(9.5.15)with respect to C.Since P(r)is a fixed polynomial,it is independent of C,hence 0=(e2+Bx+C)8%+Q()+Cx+ (9.5.18) g不市、 which can be rewritten as -Q()=(e2+Bx+C92+0+ os (9.5.19) Similarly,P()is independent of B,so differentiating(9.5.15)with respect to B gives -Q(a)=(e2+Bx+C90+0R+0s 9旦分 +B2+B (9.5.20) Now note that equation(9.5.19)matches equation(9.5.15)in form.Thus if we perform a second synthetic division of P(),ie.,a division ofQ(),yielding aremainder R+S1,then 6 OR as =-S1 (9.5.21) ac =-B1 ac To get the remaining partial derivatives,evaluate equation (9.5.20)at the two roots of the quadratic,x+and x_.Since 俗 custser g Q(x±)=R1x±+S1 (9.5.22) Numer we get 8 OR s -43126 ++那 =-x+(R1x++S1) (9.5.23) OR -+那 =-x-(R1x-+S1) (9.5.24) Solve these two equations for the partial derivatives,using E+十x-=-B T+_=C (9.5.25) and find OR OB =BR1-S1 oS aB=CR (9.5.26) Bairstow's method now consists of using Newton-Raphson in two dimensions (which is actually the subject of the next section)to find a simultaneous zero of R and S.Synthetic division is used twice per cycle to evaluate R,S and their partial derivatives with respect to B,C.Like one-dimensional Newton-Raphson,the method works well in the vicinity of a root pair(real or complex),but it can fail miserably when started at a random point.We therefore recommend it only in the context of polishing tentative complex roots
9.5 Roots of Polynomials 377 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 + b 2 ) ≡ x2 + Bx + C (9.5.14) In general if we divide a polynomial by a quadratic factor, there will be a linear remainder P(x)=(x2 + Bx + C)Q(x) + Rx + S. (9.5.15) Given B and C, R and S can be readily found, by polynomial division (§5.3). We can consider R and S to be adjustable functions of B and C, and they will be zero if the quadratic factor is a divisor of P(x). In the neighborhood of a root a first-order Taylor series expansion approximates the variation of R, S with 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) To evaluate the partial derivatives, consider the derivative of (9.5.15) with respect to C. Since P(x) is a fixed polynomial, it is independent of C, hence 0=(x2 + Bx + C) ∂Q ∂C + Q(x) + ∂R ∂C x + ∂S ∂C (9.5.18) which can be rewritten as −Q(x)=(x2 + Bx + C) ∂Q ∂C + ∂R ∂C x + ∂S ∂C (9.5.19) Similarly, P(x) is independent of B, so differentiating (9.5.15) with respect to B gives −xQ(x)=(x2 + Bx + C) ∂Q ∂B + ∂R ∂B x + ∂S ∂B (9.5.20) Now note that equation (9.5.19) matches equation (9.5.15) in form. Thus if we perform a second synthetic division of P(x), i.e., a division of Q(x), yielding a remainder 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+ and x−. Since Q(x±) = R1x± + S1 (9.5.22) we get ∂R ∂B x+ + ∂S ∂B = −x+(R1x+ + S1) (9.5.23) ∂R ∂B x− + ∂S ∂B = −x−(R1x− + S1) (9.5.24) Solve these two equations for the partial derivatives, using x+ + x− = −B x+x− = C (9.5.25) and find ∂R ∂B = BR1 − S1 ∂S ∂B = CR1 (9.5.26) Bairstow’s method now consists of using Newton-Raphson in two dimensions (which is actually the subject of the next section) to find a simultaneous zero of R and S. Synthetic division is used twice per cycle to evaluate R, S and their partial derivatives with respect to B,C. Like one-dimensional Newton-Raphson, the method works well in the vicinity of a root pair (real or complex), but it can fail miserably when started at a random point. We therefore recommend it only in the context of polishing tentative complex roots
378 Chapter 9. Root Finding and Nonlinear Sets of Equations #include #include "nrutil.h" #define ITMAX 20 At most ITMAX iterations. #define TINY 1.0e-6 void qroot(float p[],int n,float *b,float c,float eps) Given n+1 coefficients p[0..n]of a polynomial of degree n,and trial values for the coefficients of a quadratic factor x*x+b*x+c,improve the solution until the coefficients b,c change by less than eps.The routine poldiv $5.3 is used. void poldiv(float u[],int n,float v[],int nv,float g[],float r[]); int iter; float sc,sb,s,rc,rb,r,dv,delc,delb; float *q,*qq,*rem; float d[3]; granted for 19881992 q=vector(0,n); qq=vector(0,n); 1.800 rem=vector(0,n) d[2]=1.0; 872 for (iter=1;iter<=ITMAX;iter++) d[1]=(*b); d[0]=(*c); poldiv(p,n,d,2,q,rem); s=rem[0]; First division r,s. (Nort r=rem[1]; poldiv(q,(n-1),d,2,qq,rem); server computer, America to make one paper University Press. 金 THE sb=-(*c)*(rc=-rem[1]); Second division partial r,s with respect to ART rb=-(*b)*rc+(sc=-rem[0]); dv=1.0/(sb*rc-sc*rb); Solve 2x2 equation delb=(r*sc-s*rc)*dv; Programs delc=(-r*sb+s*rb)*dv; *b+=(delb=(r*sc-s*rc)*dv); OF SCIENTIFIC *c +(delc=(-r*sb+s*rb)*dv); ictly prohibited if ((fabs(delb)<=eps*fabs(*b)II fabs(*b)<TINY) &&(fabs(delc)<=eps*fabs(*c)II fabs(*c)<TINY)){ to dir free_vector(rem,0,n); Coefficients converged. free_vector(qq,0,n); free_vector(q,0,n); return; 1881892 COMPUTING (ISBN nrerror("Too many iterations in routine qroot"); v@cambri .Further reproduction, Numerical Recipes 10-621 43108 We have already remarked on the annoyance of having two tentative roots collapse to one value under polishing.You are left not knowing whether your polishing procedure has lost a root,or whether there is actually a double root, (outside 膜 which was split only by roundoff errors in your previous deflation.One solution North oftware. is deflate-and-repolish;but deflation is what we are trying to avoid at the polishing stage.An alternative is Maehly's procedure.Maehly pointed out that the derivative of the reduced polynomial P(x) P(x)≡ (9.5.27) (x-x1)…(x-) can be written as P/(x) P(x) P=c-)-女-)c-)-） ∑x-)1 (9.5.28) 1
378 Chapter 9. Root Finding and Nonlinear Sets of Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #include "nrutil.h" #define ITMAX 20 At most ITMAX iterations. #define TINY 1.0e-6 void qroot(float p[], int n, float *b, float *c, float eps) Given n+1 coefficients p[0..n] of a polynomial of degree n, and trial values for the coefficients of a quadratic factor x*x+b*x+c, improve the solution until the coefficients b,c change by less than eps. The routine poldiv §5.3 is used. { void poldiv(float u[], int n, float v[], int nv, float q[], float r[]); int iter; float sc,sb,s,rc,rb,r,dv,delc,delb; float *q,*qq,*rem; float d[3]; q=vector(0,n); qq=vector(0,n); rem=vector(0,n); d[2]=1.0; for (iter=1;iter<=ITMAX;iter++) { d[1]=(*b); d[0]=(*c); poldiv(p,n,d,2,q,rem); s=rem[0]; First division r,s. r=rem[1]; poldiv(q,(n-1),d,2,qq,rem); sb = -(*c)*(rc = -rem[1]); Second division partial r,s with respect to rb = -(*b)*rc+(sc = -rem[0]); c. dv=1.0/(sb*rc-sc*rb); Solve 2x2 equation. delb=(r*sc-s*rc)*dv; delc=(-r*sb+s*rb)*dv; *b += (delb=(r*sc-s*rc)*dv); *c += (delc=(-r*sb+s*rb)*dv); if ((fabs(delb) <= eps*fabs(*b) || fabs(*b) < TINY) && (fabs(delc) <= eps*fabs(*c) || fabs(*c) < TINY)) { free_vector(rem,0,n); Coefficients converged. free_vector(qq,0,n); free_vector(q,0,n); return; } } nrerror("Too many iterations in routine qroot"); } We have already remarked on the annoyance of having two tentative roots collapse to one value under polishing. You are left not knowing whether your polishing procedure has lost a root, or whether there is actually a double root, which was split only by roundoff errors in your previous deflation. One solution is deflate-and-repolish; but deflation is what we are trying to avoid at the polishing stage. An alternative is Maehly’s procedure. Maehly pointed out that the derivative of the reduced polynomial Pj (x) ≡ P(x) (x − x1)···(x − xj ) (9.5.27) can be written as P j (x) = P (x) (x − x1)···(x − xj ) − P(x) (x − x1)···(x − xj ) j i=1 (x − xi) −1 (9.5.28)
