c9-5 (779535), страница 2
Текст из файла (страница 2)
. . , n(9.5.8)Then we can express (9.5.6), (9.5.7) as1 n−1+=Gab1n−1+ 2 =Ha2b(9.5.9)(9.5.10)which yields as the solution for aa=G±pn(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 morerigorous justification of equation 9.5.11 is in [4].)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).Laguerre’s method is by far the most straightforward of these general, complexmethods. 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 istheoretically proved about the method’s convergence.
Much empirical experience,however, suggests that nonconvergence is extremely unusual, and, further, can almostalways be fixed by a simple scheme to break a nonconverging limit cycle. (This isimplemented in our routine, below.) An example of a polynomial that requires thiscycle-breaking scheme is one of high degree (>∼ 20), with all its roots just outside ofthe complex unit circle, approximately equally spaced around it.
When the methodconverges 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 nodisadvantage, since the polynomial itself may have complex coefficients.To motivate (although not rigorously derive) the Laguerre formulas we can notethe following relations between the polynomial and its roots and derivatives9.5 Roots of Polynomials373#include <math.h>#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 withMR different fractional values, once every MT steps, for MAXIT total allowed iterations.void laguer(fcomplex a[], int m, fcomplex *x, int *its)PiGiven the degree m and the m+1 complex coefficients a[0..m] of the polynomial mi=0 a[i]x ,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 iterationstaken 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<=MAXIT;iter++) {Loop over iterations up to allowed maximum.*its=iter;b=a[m];err=Cabs(b);d=f=Complex(0.0,0.0);abx=Cabs(*x);for (j=m-1;j>=0;j--) {Efficient computation of the polynomial andf=Cadd(Cmul(*x,f),d);its first two derivatives.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) <= err) return;We are on the root.g=Cdiv(d,b);The generic case: use Laguerre’s formula.g2=Cmul(g,g);h=Csub(g2,RCmul(2.0,Cdiv(f,b)));sq=Csqrt(RCmul((float) (m-1),Csub(RCmul((float) m,h),g2)));gp=Cadd(g,sq);gm=Csub(g,sq);abp=Cabs(gp);abm=Cabs(gm);if (abp < abm) gp=gm;dx=((FMAX(abp,abm) > 0.0 ? Cdiv(Complex((float) m,0.0),gp): RCmul(1+abx,Complex(cos((float)iter),sin((float)iter)))));x1=Csub(*x,dx);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 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 issufficiently small.The following routine implements the Laguerre method to find one root of agiven polynomial of degree m, whose coefficients can be complex. As usual, the firstcoefficient a[0] is the constant term, while a[m] is the coefficient of the highestpower of x. The routine implements a simplified version of an elegant stoppingcriterion due to Adams [5], which neatly balances the desire to achieve full machineaccuracy, on the one hand, with the danger of iterating forever in the presence ofroundoff error, on the other.374Chapter 9.Root Finding and Nonlinear Sets of Equationsif (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).}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)PiGiven 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;}}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).}nrerror("too many iterations in laguer");Very unusual — can occur only for complex roots.
Try a different starting guess for theroot.return;3759.5 Roots of PolynomialsEigenvalue Methodsam−1− am10A=...0am−2− am0101· · · − aam···0···0···10− aam00...(9.5.12)0is equivalent to the general polynomialP (x) =mXai 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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
