[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