Zeroes of a complex polynomial - Method of Jenkins and Traub

Directly to the input form

• In order to find zeroes of a complex polynomial p(z) one can consider this as the characteristic polynomial of the so called Frobenius companion matrix. This can be written differently, for example a matrix with a first subdiagonal of ones, the coefficients, divided by the leading one and taken negatively, in the last column, beginning with the lowest in row one, and zero elements otherwise: for example x³-6 x²+11 x -8 corresponds to
  

  [	0	0	8	]
  [	1	0	-11	]
  [	0	1	6	]

  The QR method for determining eigenvalues of matrices which is given by the iteration
  A_k+1=Q^T_kA_kQ_k
  where
  Q_k(A_k-s_kI)=R_k
  with Q_k unitary and R_k upper triangular, converges for well chosen shifts s_k practically always, if the matrix, as in this case, is nondecomposing upper Hessenberg. Convergence means here that the last subdiagonal element goes to zero, with an eigenvalue (a zero of p) isolated in position (n,n), such that the dimension can be reduced by one, finding the next zero in the same way.
• The method of Jenkins and Traub applies this QR technique in order to find all eigenvalues of the companion matrix. All computations necessary can be retranslated into operations on polynomials such that no matrix operations are necessary here.

Input

• The degree n of p.
• You then have the choice either to specify n complex zeroes and to compute p from these, looking whether the code delivers your input back or to specify n+1 complex coefficients a₀, ... ,a_n, with a₀ as the leading coefficient. Complex numbers must be written as pairs of real numbers in parentheses (re,im) for re + i*im, separated by comma or blank. Important: we require a_n*a₀ < > 0.
• Moreoever you may want a contourplot of |p|. In this case you must specify a border width b for a box enclosing all zeroes. |p| will be plotted over [re_min-b,re_max+b] × [im_min-b,im_max+b] where re/im_min, re/im_max are the minimal resp. maximal real resp. imaginary part of a zero. Avoid a large b since outside the region of zeroes the value of p grows very rapidly and you will see little useful.

Output

• The computed zeroes.
• If desired the contourplot of |p|. The area will be adapted internally if necessary, for example if all zeroes are real.

Questions?!

• What can you learn about sensitivity of zeroes with respect to coefficient perturbations, depending on the distribution of the zeroes (linear or geometric, on the unit circle) and what can you observe when there are multiple zeroes, for example (x-1)⁴ and (x-1)⁴+(x²-2x+0.99)*eps?
• Does it indeed ''never fail'' ?
• Did you believe that a ''harmless function'' like a polynomial would look like this?

To the input form