numerical stability of root identification via Newton-Raphson iteration of Stieltjes residue sums

I have asked several questions on math.SE in order to compute numerically the poles of high-degree Padé approximations for ex, because a computation directly from the polynomial coefficients has poor numerical stability. I finally blundered upon an approach alluded to a paper by Campos and Calderón and a related item mentioned in a paper by Grünbaum which states that for polynomials with simple roots

If
P(x)=(xx1)(xx2)(xxn)
then at any root xk of P we have
P

In other words,

\frac{P”(x_k)}{2P'(x_k)} = \sum\limits_{j \neq k}\frac{1}{x_k-x_j}

(I called this a “Stieltjes residue sum” in the title since it seems like his name pops up in this arena, and Grünbaum’s article cites some work he did in electrostatics applications with these sorts of sums. If there’s a more accurate name, please correct me + I’ll change the title.)

For certain polynomials P(x) satisfying a nonlinear second-order differential equation, the left side becomes very simple. For example, for the Bessel polynomials y_n(x) and \theta_n(x) with roots y_k and z_k, respectively,

\frac{1}{y_j}+\frac{1}{y_j^2} = \sum\limits_{k\neq j} \frac{1}{y_k – y_j}

-1-\frac{n}{z_j} = \sum\limits_{k\neq j} \frac{1}{z_k – z_j}

and for the denominator D_{pq}(x) = {}_1F_1(-q;-p-q;x) of the Padé approximation of e^{-x} with roots x_k,

\frac{x_k+p+q}{2x_k} = \sum\limits_{k\neq j}\frac{1}{x_k-x_j}

These can be turned into nonlinear vector equations \mathbf F(\mathbf x) = \mathbf 0 where \mathbf x is a vector of roots, and solved iteratively by the Newton-Raphson method given an initial guess that is vaguely close to the exact values. I’ve done this for the Bessel polynomial \theta_n(x) (described in detail on math.stackexchange.com) and for D_{pq}(x) and got very good results with fast convergence (8-10 iterations) even for 500th-degree polynomials.

What I don’t understand is why this works so well. Can anyone explain why this form of a polynomial would converge so well? Finding roots of the polynomials given the coefficients is a very ill-conditioned problem. And this method has n! possible solutions, so it seems like any iterative approach would be fraught with instability.

Answer

Attribution
Source : Link , Question Author : Jason S , Answer Author : Community

Leave a Comment