In [1]:
from scipy import *;
from scipy.linalg import eigh
In [2]:
k =mat([[ 54, -34, 8, 0, 0],
[ -34, 48, -30, 7, 0],
[ 8, -30, 42, -26, 6],
[ 0, 7, -26, 31, -12],
[ 0, 0, 6, -12, 6]])*0.65
mass = mat(eye(5)) ; mass[4,4] = 4 ; mass = mass*1.0
In [3]:
evals, evecs = eigh(k, mass)
evecs = mat(evecs)
w2 = mat(diag(evals)) ; w1 = mat(diag(sqrt(evals)))
print evals
print sqrt(evals)
print 1/(1-1/evals)
print evecs
In [4]:
load = mat('0;0;0;2;-1')
qload = evecs.T*load ; Gamma = diag(asarray(qload).flatten())
print Gamma
In [5]:
R = mass*evecs*Gamma
print R
In [6]:
print [(R[i,:]).sum() for i in range(5)]
print [(R[:,i]).sum() for i in range(5)]
In [7]:
set_printoptions(precision=6, linewidth=120, suppress=1)
f = copy(k.I)
print "Now, the modal flexibility matrices, F_i = \psi*\psi^T/\omega^2,\n"
for i in range(5):
print "F_%d ="%(i+1,)
f_i = evecs[:,i]*evecs[:,i].T/evals[i]
print f_i
f = f - f_i
print
print 'The following matrix, F-sum(F_i), should be close to zero...'
In [9]:
print f
print "... and so it seems."
In [10]:
print "The correction matrix for modes 3,4,5"
f = k.I - evecs[:,0]*evecs[:,0].T/evals[0] - evecs[:,1]*evecs[:,1].T/evals[1]
print f,"\n"
print "x = \psi_1 q_1 + \psi_2 q_2 + ((F-F_1-F_2)*r) f(t)"
print " ^^^^^^^^^^^^^^^"
print " constant"
let's try $p(t) = \boldsymbol{r} \\,\sin(t)$, that is $\omega=1.0$
In [ ]:
In [ ]: