We take a string, divide it up into N points, and then approximate Hooke's law between the points. The eigenvalues of the resulting system of equations are the modes of vibration
# TODO work this out better!
# NOTE there might be an off-by-sign error here
N = 25
A = zeros(N,N)
A = A + diagm(-2*ones(N))
A = A + diagm(1*ones(N-1),1)
A = A + diagm(1*ones(N-1),-1)
using Plots
# We actually want the smallest eigenvalues
# So we are going to multiply A by -1
# to make the order that Julia returns them
# correct
vals,vecs = eig(-A)
plot(vecs[:,1:6],layout=6)
## Create a symmetric matrix with largest eigenvalue 1
A = [2.8 0.6; 0.6 1.2]/3
println("Eigenvalues")
@show eigvals(A)
println("Eigenvectors = ")
println( eigvecs(A) )
## Show that powers of the matrix converge to the eigenvector
println("A^25 = ")
println(A^25)
; # hide last output
vecs = eigvecs(A)
vecs[:,2]*vecs[:,2]'
y = A^25*randn(2,1)
y/norm(y)
## Try out a practical version of the power method
x = randn(2,1)
for i=1:42
y = A*x
x = y/norm(y)
@printf(" x[1] = %.3f x[2] = %.3f \n", x[1], x[2])
end
## Try out the same method on the matrix we first had
A = [1 -10; -5 -4.]
x = randn(2,1)
for i=1:100
y = A*x
x = y/norm(y)
@printf(" x[1] = %.3f x[2] = %.3f \n", x[1], x[2])
end