c9-5 (779535), страница 3
Текст из файла (страница 3)
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[])PmiFind 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];Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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 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 matrix376Chapter 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;}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 advantageSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).}free_matrix(hess,1,MAXM,1,MAXM);3779.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)R(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)we get∂R∂Sx− += −x− (R1 x− + S1 )∂B∂BSolve these two equations for the partial derivatives, usingx+ + x− = −Bx+ x− = C(9.5.24)(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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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, C378Chapter 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.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 asXP (x)P 0 (x)−(x − xi )−1 (9.5.28)(x − x1 ) · · · (x − xj ) (x − x1 ) · · · (x − xj )jPj0 (x) =i=1Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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];3799.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 0 (xP (xk )Pj−1k ) − P (xk )i=1 (xk − xi )(9.5.29)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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
