Thought I'd do a self-plug and mention that you can do this and other
polynomial manipulations in ApproxFun.jl (which uses a colleague matrix and
gets 10 digits accuracy for this example):
using ApproxFun
f=Fun(x->prod([1.0:20.0]-x),[0,20],21)
roots(f)
On Tuesday, May 6, 2014 at 6:14:44 AM UTC+10, Hans W Borchers wrote:
>
> The following is the so-called Wilkinson polynomial, that has zeros at 1,
> 2, ..., 20 :
>
> julia> using Polynomial
>
>
> julia> pW = poly([1.0:20.0])
>
>
> julia> roots(pW)
> 20-element Array{Complex{Float64},1}:
> 1.0-0.0im
> 2.0-0.0im
> 3.0-0.0im
> 3.99999-0.0im
> 5.00006-0.0im
> 6.00181-0.0im
> 6.95031-0.0im
> 8.09172-0.579086im
> 8.09172+0.579086im
> 9.70891-1.64679im
> 9.70891+1.64679im
> 11.7999-2.59604im
> 11.7999+2.59604im
> 14.3615-3.12883im
> 14.3615+3.12883im
> 17.1543-2.89387im
> 17.1543+2.89387im
> 19.5758-1.70785im
> 19.5758+1.70785im
> 20.6635-0.0im
>
>
> Compare this with the same computation in Matlab:
>
> >> pw = poly(1:20);
> >> roots(pw)
>
> ans =
>
> 20.0003
> 18.9972
> 18.0112
> 16.9711
> 16.0483
> 14.9354
> 14.0653
> 12.9491
> 12.0334
> 10.9840
> 10.0061
> 8.9984
> 8.0003
> 7.0000
> 6.0000
> 5.0000
> 4.0000
> 3.0000
> 2.0000
> 1.0000
>
>
> The roots of this polynomial are known to be difficult to determine, but
> the results in Julia's Polynomial package seem really far off.
>
