Jenkins–Traub algorithm: Difference between revisions

Revision as of 23:41, 12 December 2010

Revision as of 10:04, 29 December 2010

Line 1:

The '''Jenkins-Traub algorithm for polynomial zeros''' is a fast globally convergent iterative method. It has been described as ''practically a standard in black-box polynomial root-finders''.<ref>Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (2007), Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, page 470.</ref>.

Given a polynomial ''P'',

Line 9:

==Overview==

The Jenkins-Traub algorithm calculates all of the roots of a [[polynomial]] with complex coefficients. The algorithm starts by checking the polynomial for the occurrence of very large or very small roots. If necessary, the coefficients are rescaled by a rescaling of the variable. In the algorithm proper, roots are found one by one and generally in increasing size. After each root is found, the polynomial is deflated by dividing off the corresponding linear factor. Indeed, the factorization of the polynomial into the linear factor and the remaining deflated polynomial is already a result of the root-finding procedure. The root-finding procedure has three stages that correspond to different variants of the [[inverse power iteration]]. See Jenkins and [[Joseph F Traub|Traub]].<ref>Jenkins, M. A. and Traub, J. F. (1970), [http://www.springerlink.com/content/q6w17w30035r2152/?p=ae17d723839045be82d270b45363625f&pi=1 A Three-Stage Variables-Shift Iteration for Polynomial Zeros and Its Relation to Generalized Rayleigh Iteration], Numer. Math. 14, 252-263.</ref> .

A description can also be found in Ralston and [[Philip Rabinowitz (mathematician)|

Rabinowitz]]<ref>Ralston, A. and Rabinowitz, P. (1978), A First Course in Numerical Analysis, 2nd ed., McGraw-Hill, New York.</ref> p.&nbsp;383.

The algorithm is similar in spirit to the two-stage algorithm studied by Traub.<ref>Traub, J. F. (1966), [http://links.jstor.org/sici?sici=0025-5718(196601)20%3A93%3C113%3AACOGCI%3E2.0.CO%3B2-3 A Class of Globally Convergent Iteration Functions for the Solution of Polynomial Equations], Math. Comp., 20(93), 113-138.</ref>.

=== Root-finding procedure ===

Line 66:

Now choose <math>s=R\cdot \exp(i\,\phi_\text{random})</math> on the circle of this radius. The sequence of polynomials <math>H^{(\lambda+1)}(z)</math>, <math>\lambda=M,M+1,\dots,L-1</math>, is generated with the fixed shift value <math>s_\lambda=s</math>. During this iteration, the current approximation for the root

:<math>t_\lambda=s-\frac{P(s)}{\bar H^{(\lambda)}(s)}</math>

is traced. The second stage is finished successfully if the conditions

:<math>

Line 72:

</math> and <math>

|t_\lambda-t_{\lambda-1}|<\tfrac12\,|t_{\lambda-1}|

</math>

are simultaneously met. If there was no success after some number of iterations, a different random point on the circle is tried. Typically one uses a number of 9 iterations for polynomials of moderate degree, with a doubling strategy for the case of multiple failures.

Line 104:

is precisely a Newton-Raphson iteration performed on certain [[rational functions]]. More precisely, Newton-Raphson is being performed on a sequence of rational functions

:<math>W^\lambda(z)=\frac{P(z)}{H^\lambda(z)}</math>.

For <math>\lambda</math> sufficiently large,

Line 112:

is as close as desired to a first degree polynomial

:<math>z-\alpha_1</math>,

where <math>\alpha_1</math> is one of the zeros of <math>P</math>. Even though Stage 3 is precisely a Newton-Raphson iteration, differentiation is not performed.

Line 172:

\left|\frac{\alpha_1-s}{\alpha_2-s}\right|^{\lambda-M}\right)

</math>

*:and

*:<math>

s-\frac{P(s)}{\bar H^{(\lambda)}(s)}

Line 184:

\right)

</math>

*:and

*:<math>

s_{\lambda+1}=

Line 200:

:<math>P(X)=(X-\alpha_1)\cdot P_1(X)</math>

with a root <math>\alpha_1\in\C</math> and <math>P_1(X)=P(X)/(X-\alpha_1)</math> the remaining factor of degree ''n-1'' as the eigenvector equation for the multiplication with the variable ''X'', followed by remainder computation with divisor ''P(X)'',

:<math>M_X(H)=(X\cdot H(X)) \bmod P(X)\,.</math>

This maps polynomials of degree at most ''n-1'' to polynomials of degree at most ''n-1''. The eigenvalues of this map are the roots of ''P(X)'', since the eigenvector equation reads

:<math>0=(M_X-\alpha\cdot id)(H)=((X-\alpha)\cdot H) \bmod P\,,</math>

Line 238:

*[http://www.netlib.org/toms/493 RPOLY] Algorithm 493 from ACM TOMS) on Netlib. In Fortran.

{{DEFAULTSORT:Jenkins-Traub Algorithm}}

[[Category:Numerical analysis]]

[[Category:Root-finding algorithms]]