[Numpy-discussion] roots and high-order polynomial
Pauli Virtanen
pav@iki...
Fri Jul 3 14:08:44 CDT 2009
On 2009-07-03, Charles R Harris <charlesr.harris@gmail.com> wrote:
> roots? The connection between polynomial coefficients and polynomial values
> becomes somewhat vague when the polynomial degree becomes large, it is
> numerically ill conditioned.
In addition to switching to higher precision than machine
precision, another alternative could be choosing a stabler (eg.
Chebyshev) polynomial basis than the {1,x,x^2,...} one.
For the Chebyshev series f(x) = sum_{j=0}^N a_j T_j(x), the
companion matrix is (check!)
A = zeros((N, N))
A[0,1] = 1
j = arange(1, N-1)
A[j, j+1] = 0.5
A[j, j-1] = 0.5
A[-1,:] = -0.5*a[:-1]/a[-1]
A[-1,-2] += 0.5
and the zeros are
x = linalg.eig(A)[0]
See eg. http://dx.doi.org/10.1016/j.apnum.2005.09.007 for more.
--
Pauli Virtanen
More information about the NumPy-Discussion
mailing list