## 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