In [1]:
macro myevalpoly(z,a...)
isempty(a) && error("You forgot to pass coefficients!")
ex = :($(a[length(a)]))
for i in 1:length(a)-1
ex = :($ex * $(z) + $(a[length(a)-i]) )
end
println(ex)
ex
end
Out[1]:
In [2]:
@myevalpoly 7 2 3 4 5
Out[2]:
In [3]:
@evalpoly 7 2 3 4 5
Out[3]:
First, we need to construct coefficients $a_k$. For the polynomial $\prod_{i=1}^4 (x-z_i)$, we have the coefficients $$\left(
\begin{array}{c}
z_1 z_2 z_3 z_4 \\
-z_1 z_2 z_3-z_1 z_4 z_3-z_2 z_4 z_3-z_1 z_2 z_4 \\
z_1 z_2+z_3 z_2+z_4 z_2+z_1 z_3+z_1 z_4+z_3 z_4 \\
-z_1-z_2-z_3-z_4 \\
1 \\
\end{array}
\right),$$ thus we can exploit the structure and write a double for loop to calculate the coefficients. A more general formula is
$$
\begin{cases}
1 = a_{n}\\
x_{1}+x_{2}+\dots +x_{n-1}+x_{n}=-a_{n-1}\\
(x_{1}x_{2}+x_{1}x_{3}+\cdots +x_{1}x_{n})+(x_{2}x_{3}+x_{2}x_{4}+\cdots +x_{2}x_{n})+\cdots +x_{n-1}x_{n}=a_{n-2}\\
\quad \vdots \\
x_{1}x_{2}\dots x_{n}=(-1)^{n}a_{0}.\end{cases}
$$
Checkout Vieta's formulas for more information.
In [4]:
function root2coeff(z::AbstractVector{T}) where T
N = length(z)
co = zeros(T, N+1)
# The last coefficient is always one
co[end] = 1
# The outer loop adds one root at a time
for j in 1:N, i in j:-1:1
co[end-i] -= z[j]*co[end-i+1]
end
co
end
@show typemax(Int), typemax(Int128)
root2coeff(1:20)
Out[4]:
Those numbers are close to typemax(Int), so integer overflows may occur, lets use Int128 instead.
In [5]:
root2coeff(Int128(1):20)
Out[5]:
Next, we need to construct a companion matrix and solve for roots. A companion matrix is in the form of $$ \begin{bmatrix}0&0&\dots &0&-z_{1}\\1&0&\dots &0&-z_{2}\\0&1&\dots &0&-z_{3}\\\vdots &\vdots &\ddots &\vdots &\vdots \\0&0&\dots &1&-z_{n-1}\end{bmatrix}. $$
In [6]:
function poly_roots(z)
len = length(z)
# construct the ones part
mat = diagm(ones(len-2), -1)
# insert coefficients
mat[:, end] = -z[1:end-1]
eigvals(mat)
end
Out[6]:
We have everything ready now. We just need to calculate all the roots and plot it.
In [7]:
using Random
Random.seed!(1)
function wilkinson_poly_roots(n=100)
# original coefficients
coeff = root2coeff(Int128(1):20)
rts = Vector{Complex{Float64}}[]
# add perturbation
for i in 1:n
pert_coeff = coeff.*(1+rand(21)*1e-10)
push!(rts, poly_roots(pert_coeff))
end
rts
end
using Plots; gr()
function plt_wilkinson_roots(rts)
# plot roots without perturbation
plt = scatter(1:20, zeros(20), color = :green, markersize = 5, legend=false)
for i in eachindex(rts)
# plot roots with perturbation
scatter!(plt, real.(rts[i]), imag.(rts[i]), color = :red, markersize = .5)
end
plt
end
wilkinson_poly_roots() |> plt_wilkinson_roots