Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 89
Текст из файла (страница 89)
Try a different starting guess for theroot.return;}Here is a driver routine that calls laguer in succession for each root, performsthe deflation, optionally polishes the roots by the same Laguerre method — if youare not going to polish in some other way — and finally sorts the roots by their realparts. (We will use this routine in Chapter 13.)#include <math.h>#include "complex.h"#define EPS 2.0e-6#define MAXM 100A small number, and maximum anticipated value of m.void zroots(fcomplex a[], int m, fcomplex roots[], int polish)"iGiven the degree m and the m+1 complex coefficients a[0..m] of the polynomial mi=0 a(i)x ,this routine successively calls laguer and finds all m complex roots in roots[1..m].
Theboolean 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<=m;j++) ad[j]=a[j];Copy of coefficients for successive deflation.for (j=m;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) <= 2.0*EPS*fabs(x.r)) x.i=0.0;roots[j]=x;b=ad[j];Forward deflation.for (jj=j-1;jj>=0;jj--) {c=ad[jj];ad[jj]=b;b=Cadd(Cmul(x,b),c);}}if (polish)for (j=1;j<=m;j++)Polish the roots using the undeflated coeffilaguer(a,m,&roots[j],&its);cients.for (j=2;j<=m;j++) {Sort roots by their real parts by straight inx=roots[j];sertion.for (i=j-1;i>=1;i--) {if (roots[i].r <= x.r) break;roots[i+1]=roots[i];}roots[i+1]=x;}}3759.5 Roots of PolynomialsEigenvalue MethodsThe 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 notgenerally an efficient way to find eigenvalues. Turning matters around, we canuse the more efficient eigenvalue methods that are discussed in Chapter 11 to findthe roots of arbitrary polynomials. You can easily verify (see, e.g., [6]) that thecharacteristic polynomial of the special m × m companion matrixam−1− am10A=...0am−2− am0101· · · − aam···0···0···10− aam00...(9.5.12)0is equivalent to the general polynomialP (x) =mai x i(9.5.13)i=0If the coefficients ai are real, rather than complex, then the eigenvalues of A can befound using the routines balanc and hqr in §§11.5–11.6 (see discussion there).
Thismethod, implemented in the routine zrhqr following, is typically about a factor 2slower than zroots (above). However, for some classes of polynomials, it is a morerobust technique, largely because of the fairly sophisticated convergence methodsembodied in hqr. If your polynomial has real coefficients, and you are havingtrouble with zroots, then zrhqr is a recommended alternative.#include "nrutil.h"#define MAXM 50void zrhqr(float a[], int m, float rtr[], float rti[])"miFind all the roots of a polynomial with real coefficients,i=0 a(i)x , given the degree mand the coefficients a[0..m]. The method is to construct an upper Hessenberg matrix whoseeigenvalues are the desired roots, and then use the routines balanc and hqr. The real andimaginary 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];376Chapter 9.Root Finding and Nonlinear Sets of Equationsfor (k=j-1;k>=1;k--) {if (rtr[k] <= xr) break;rtr[k+1]=rtr[k];rti[k+1]=rti[k];}rtr[k+1]=xr;rti[k+1]=xi;}free_matrix(hess,1,MAXM,1,MAXM);}Other Sure-Fire TechniquesThe Jenkins-Traub method has become practically a standard in black-boxpolynomial root-finders, e.g., in the IMSL library [3].
The method is too complicatedto discuss here, but is detailed, with references to the primary literature, in [4].The Lehmer-Schur algorithm is one of a class of methods that isolate roots inthe complex plane by generalizing the notion of one-dimensional bracketing. It ispossible to determine efficiently whether there are any polynomial roots within acircle of given center and radius. From then on it is a matter of bookkeeping tohunt down all the roots by a series of decisions regarding where to place new trialcircles. Consult [1] for an introduction.Techniques for Root-PolishingNewton-Raphson works very well for real roots once the neighborhood ofa root has been identified. The polynomial and its derivative can be efficientlysimultaneously evaluated as in §5.3.
For a polynomial of degree n with coefficientsc[0]...c[n], the following segment of code embodies one cycle of NewtonRaphson:p=c[n]*x+c[n-1];p1=c[n];for(i=n-2;i>=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 thecomplex roots, either directly, or by looking for quadratic factors.Direct polishing by Newton-Raphson is straightforward for complex roots if theabove code is converted to complex data types. With real polynomial coefficients,note that your starting guess (tentative root) must be off the real axis, otherwiseyou will never get off that axis — and may get shot off to infinity by a minimumor maximum of the polynomial.For real polynomials, the alternative means of polishing complex roots (or, for thatmatter, double real roots) is Bairstow’s method, which seeks quadratic factors.
The advantage3779.5 Roots of Polynomialsof going after quadratic factors is that it avoids all complex arithmetic. Bairstow’s methodseeks a quadratic factor that embodies the two roots x = a ± ib, namelyx2 − 2ax + (a2 + b2 ) ≡ x2 + Bx + C(9.5.14)In general if we divide a polynomial by a quadratic factor, there will be a linear remainderP (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 canconsider R and S to be adjustable functions of B and C, and they will be zero if thequadratic factor is zero.In the neighborhood of a root a first-order Taylor series expansion approximates thevariation of R, S with respect to small changes in B, CR(B + δB, C + δC) ≈ R(B, C) +∂R∂RδB +δC∂B∂C(9.5.16)∂S∂SδB +δC(9.5.17)∂B∂CTo evaluate the partial derivatives, consider the derivative of (9.5.15) with respect to C.
SinceP (x) is a fixed polynomial, it is independent of C, henceS(B + δB, C + δC) ≈ S(B, C) +0 = (x2 + Bx + C)∂Q∂R∂S+ Q(x) +x+∂C∂C∂C(9.5.18)which can be rewritten as∂Q∂R∂S+x+(9.5.19)∂C∂C∂CSimilarly, P (x) is independent of B, so differentiating (9.5.15) with respect to B gives−Q(x) = (x2 + Bx + C)∂R∂S∂Q+x+(9.5.20)∂B∂B∂BNow note that equation (9.5.19) matches equation (9.5.15) in form. Thus if we perform asecond synthetic division of P (x), i.e., a division of Q(x), yielding a remainder R1 x+S1 , then−xQ(x) = (x2 + Bx + C)∂R∂S(9.5.21)= −R1= −S1∂C∂CTo get the remaining partial derivatives, evaluate equation (9.5.20) at the two roots of thequadratic, x+ and x− .
SinceQ(x± ) = R1 x± + S1(9.5.22)∂R∂Sx+ += −x+ (R1 x+ + S1 )∂B∂B(9.5.23)∂R∂Sx− += −x− (R1 x− + S1 )∂B∂BSolve these two equations for the partial derivatives, using(9.5.24)we getx+ + x− = −Bx+ x− = C(9.5.25)and find∂S∂R(9.5.26)= BR1 − S1= CR1∂B∂BBairstow’s method now consists of using Newton-Raphson in two dimensions (which isactually the subject of the next section) to find a simultaneous zero of R and S. Syntheticdivision is used twice per cycle to evaluate R, S and their partial derivatives with respect toB, C. Like one-dimensional Newton-Raphson, the method works well in the vicinity of a rootpair (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.378Chapter 9.Root Finding and Nonlinear Sets of Equations#include <math.h>#include "nrutil.h"#define ITMAX 20#define TINY 1.0e-6At most ITMAX iterations.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 coefficientsof a quadratic factor x*x+b*x+c, improve the solution until the coefficients b,c change by lessthan 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 torb = -(*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 rootscollapse to one value under polishing.
You are left not knowing whether yourpolishing 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 solutionis deflate-and-repolish; but deflation is what we are trying to avoid at the polishingstage. An alternative is Maehly’s procedure. Maehly pointed out that the derivativeof the reduced polynomialPj (x) ≡P (x)(x − x1 ) · · · (x − xj )(9.5.27)can be written asPj (x)jP (x)P (x)−=(x − xi )−1 (9.5.28)(x − x1 ) · · · (x − xj ) (x − x1 ) · · · (x − xj )i=13799.6 Newton-Raphson Method for Nonlinear Systems of EquationsHence one step of Newton-Raphson, taking a guess xk into a new guess xk+1 ,can be written asxk+1 = xk −P (xP (xk )"j−1k ) − P (xk )i=1 (xk − xi )(9.5.29)This equation, if used with i ranging over the roots already polished, will prevent atentative root from spuriously hopping to another one’s true root.
It is an exampleof so-called zero suppression as an alternative to true deflation.Muller’s method, which was described above, can also be useful at thepolishing stage.CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 7. [1]Peters G., and Wilkinson, J.H. 1971, Journal of the Institute of Mathematics and its Applications,vol. 8, pp.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
