Vibration analysis of a string with eigenvalues

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

In [21]:
# 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)
Out[21]:
25x25 Array{Float64,2}:
 -2.0   1.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  1.0  -2.0   1.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   1.0  -2.0   1.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   1.0  -2.0   1.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   1.0  -2.0   1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   1.0  -2.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      1.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0     -2.0   1.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   1.0  -2.0   1.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   1.0  -2.0   1.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   1.0  -2.0   1.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   1.0  -2.0   1.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   1.0  -2.0
In [19]:
using Plots
In [22]:
# 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)
[Plots.jl] Initializing backend: pyplot
WARNING: No working GUI backend found for matplotlib.
Out[22]:
In [ ]:

The power method

In [23]:
## 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
Eigenvalues
eigvals(A) = [0.3333333333333333,0.9999999999999999]
Eigenvectors = 
[0.31622776601683794 -0.9486832980505138
 -0.9486832980505138 -0.31622776601683794]
A^25 = 
[0.9000000000001165 0.29999999999964533
 0.2999999999996454 0.100000000001062]
In [27]:
vecs = eigvecs(A)
vecs[:,2]*vecs[:,2]'
Out[27]:
2x2 Array{Float64,2}:
 0.9  0.3
 0.3  0.1
In [51]:
y = A^25*randn(2,1)
y/norm(y)
Out[51]:
2x1 Array{Float64,2}:
 -0.948683
 -0.316228
In [52]:
## 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
 x[1] = 0.142   x[2] = 0.990  
 x[1] = 0.615   x[2] = 0.789  
 x[1] = 0.858   x[2] = 0.514  
 x[1] = 0.923   x[2] = 0.385  
 x[1] = 0.941   x[2] = 0.340  
 x[1] = 0.946   x[2] = 0.324  
 x[1] = 0.948   x[2] = 0.319  
 x[1] = 0.948   x[2] = 0.317  
 x[1] = 0.949   x[2] = 0.317  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
 x[1] = 0.949   x[2] = 0.316  
In [18]:
## 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
 x[1] = -0.967   x[2] = -0.254  
 x[1] = 0.260   x[2] = 0.966  
 x[1] = -0.876   x[2] = -0.481  
 x[1] = 0.530   x[2] = 0.848  
 x[1] = -0.796   x[2] = -0.605  
 x[1] = 0.634   x[2] = 0.773  
 x[1] = -0.750   x[2] = -0.662  
 x[1] = 0.676   x[2] = 0.737  
 x[1] = -0.727   x[2] = -0.687  
 x[1] = 0.694   x[2] = 0.720  
 x[1] = -0.716   x[2] = -0.698  
 x[1] = 0.701   x[2] = 0.713  
 x[1] = -0.711   x[2] = -0.703  
 x[1] = 0.704   x[2] = 0.710  
 x[1] = -0.709   x[2] = -0.705  
 x[1] = 0.706   x[2] = 0.708  
 x[1] = -0.708   x[2] = -0.706  
 x[1] = 0.707   x[2] = 0.708  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707  
 x[1] = -0.707   x[2] = -0.707  
 x[1] = 0.707   x[2] = 0.707