Good catch. This looks almost entirely due to the order in which
Polynomial.jl is creating the companion matrix. The modified version
at https://github.com/loladiro/Polynomials.jl which stores the coefficients
in the mathematical order instead of the Matlab order does a little better,
but worryingly I'm getting noticeably different results between 32 bit and
64 bit.
Just transposing the companion matrix gives much better-behaved results for
this polynomial, not exactly the same as Matlab but pretty close.
On Monday, May 5, 2014 1:14:44 PM UTC-7, 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.
>
