using LinearAlgebra, PyPlot, Interact

A = randn(5,5)
A = A * Diagonal(1:5) / A

eigvals(A)

λ, X = eigen(A, sortby=λ->-abs(λ)) # sort eigenvals by -|λ|
λ

x = [5,4,3,2,1] # arbitrary initial vector
# @manipulate
for n = 1:40
    y = x
    for i = 1:n
        y = A*y
        y = y / norm(y)
    end
    display(
        hbox(vbox(HTML("<b>power iteration $n</b>"),y),
         vline(), 
         vbox(HTML("<b>eigenvector</b>"),X[:,1]))
        )
end

d = Float64[]
y = x
for i = 1:100
    y = A*y
    y = y / norm(y)
    push!(d, min(norm(y - X[:,1]), norm(-y - X[:,1]))) # pick the better of the two signs
end
semilogy(1:length(d), d, "b.-")
xlabel("number of power steps")
ylabel("error in eigenvector")
title(L"convergence of power method for $\lambda=1,2,3,4,5$")

semilogy(1:length(d), d, "b.-")
semilogy(1:length(d), (4/5).^(1:length(d)), "k--")
xlabel("number of power steps")
ylabel("error in eigenvector")
title(L"convergence of power method for $\lambda=1,2,3,4,5$")
legend(["error", L"(4/5)^n"])

y = x
for i = 1:30
    y = A*y
    y = y / norm(y)
end
y

A*y

(A*y) ./ y   # divide each element of Ay elementwise (./ in Julia) by the corresponding element of y

(y'A*y) / (y'y) # a Rayleigh quotient in Julia

Δλ = Float64[]
y = x
for i = 1:100
    y = A*y
    y = y / norm(y)
    λ̃ = (y'A*y) / (y'y)
    push!(Δλ, abs(λ̃ - 5))
end
semilogy(1:length(Δλ), Δλ, "b.-")
semilogy(1:length(Δλ), (4/5).^(1:length(Δλ)), "k--")
xlabel("number of power steps")
ylabel("error in eigenvalue")
title(L"convergence of power method $\lambda$ for $\lambda=1,2,3,4,5$")
legend(["error", L"(4/5)^n"])

Δλ = Float64[]
y = x
for i = 1:15
    y = (A - 2.1I) \ y
    y = y / norm(y)
    λ̃ = (y'A*y) / (y'y)
    push!(Δλ, abs(λ̃ - 2))
end
semilogy(1:length(Δλ), Δλ, "b.-")
semilogy(1:length(Δλ), (1/9).^(1:length(Δλ)), "k--")
xlabel("number of power steps")
ylabel("error in eigenvalue")
title(L"convergence of shift-and-invert method $\lambda$ for $\lambda=1,2,3,4,5$, $\mu = 2.1$")
legend(["error", L"1/9^n"])