I need to obtain (numerically) the roots of the denominator in the Padé approximation to the exponential function $e^{-x}$, in Python. I can calculate its coefficients in closed form (see below). But the numpy.roots() function to find roots of a polynomial of known coefficients has a little footnote:
The algorithm relies on computing the eigenvalues of the companion matrix [R279].
and the companion matrix becomes very ill-conditioned as the degree increases. (2-norm condition number of approx 2.23e+10 for a ratio of 9th-degree polynomials)
So it seems like I would be better off with a more direct approach to getting the roots; for matrix exponential $e^A$ the Padé coefficients for numerator $N_{pq}(A)$ of degree $p$ and denominator $D_{pq}(A)$ of degree $q$ for and are given in p.9 of Moler and Van Loan's paper on matrix exponentials; if I substitute $A=-x$ then I get this for the denominator:
$$D_{pq}(x) = \sum\limits_{j=0}^q\frac{(p+q-j)!q!}{(p+q)!j!(q-j)!}x^j$$
Is there a way to find the roots more directly? They are very nicely structured, I just don't know what that structure is; here is a complex-plane plot of the poles ("x") and zeros ("o") for $p=q=9$:
If I use numpy.roots() (which, again, relies on computing the eigenvalues of the companion matrix), then things seem to be ok for a while but then maybe break down...
$p=q=16$:
$p=q=27$:
$p=q=28$:
$p=q=32$:
and soon it looks like a marching band gone crazy (though symmetrically crazy in conjugate pairs).