Yes, Polynomial is using a different convention than Matlab or what you used below in how it constructs the companion matrix. Polynomials.jl uses yet another convention. Both produce more accurate (comparable to Matlab) results for the roots of the Wilkinson polynomial if you just switch the indices during construction of the companion matrix. See https://github.com/vtjnash/Polynomial.jl/blob/master/src/Polynomial.jl#L350-L353 for Polynomial, or https://github.com/loladiro/Polynomials.jl/blob/master/src/Polynomials.jl#L324-L325 for Polynomials.
I think the intent (see https://github.com/vtjnash/Polynomial.jl/issues/5) is to deprecate Polynomial and switch the coefficient order by developing under the Polynomials name going forward, but Keno's probably been too busy to register the new package, turn on issues, etc. I'm doing some work on piecewise stuff in a branch of Polynomials, I might just adopt the package. I know David de Laat has put together packages for sparse multivariate polynomials https://github.com/daviddelaat/MultiPoly.jl and orthogonal polynomials https://github.com/daviddelaat/Orthopolys.jl, don't think they're registered but it might make sense to eventually unify all of these into the same package. On Tuesday, May 6, 2014 8:16:08 PM UTC-7, Hans W Borchers wrote: > > So where is the problem/bug? > I did construct the companion matrix manually, starting from the vector of > coefficients of the Wilkinson polynomial: > > julia> pW = [ > 1.0, -210, 20615, > -1256850, 53327946, -1672280820, > 40171771630, -756111184500, 11310276995381, > -135585182899530, 1307535010540395, -10142299865511450, > 63030812099294896, -311333643161390656, 1206647803780373248, > -3599979517947607040, 8037811822645052416, -12870931245150988288, > 13803759753640704000, -8752948036761600000, 2432902008176640000]; > > > julia> pw = pW[2:21] / pW[1]; # normalize 1st > coefficient > > julia> A = vcat(-pw', hcat(eye(19), zeros(19)); # companion matrix > > julia> eigvals(A) > 20-element Array{Float64,1}: > 20.0003 > 18.9979 > 18.0072 > 16.9854 > 16.0183 > 14.9856 > 14.0067 > 12.9987 > 12.0002 > 10.9992 > 10.0009 > 8.99944 > 8.0002 > 6.99995 > 6.00001 > 5.0 > 4.0 > 3.0 > 2.0 > 1.0 > > This is almost the same result as with Matlab above. Does this mean, > Polynomial is constructing a different / wrong companion matrix? > > > On Monday, May 5, 2014 11:24:47 PM UTC+2, Tony Kelman wrote: >> >> 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. >> >> >>
