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




[  2.21379558e-02   1.47531594e+00   1.14460116e+01   3.52889322e+01
6.64926023e+01]
[ 0.14878829  1.21462584  3.38319548  5.94044882  8.15429962]
[-0.02263914  3.10386382  1.09573032  1.02916393  1.0152689 ]
[[ 0.03861265 -0.26504416 -0.52623003  0.60662458 -0.53230558]
[ 0.11061399 -0.54887627 -0.51096499 -0.12909326  0.63933488]
[ 0.20963847 -0.61887301  0.22213708 -0.52660579 -0.49637624]
[ 0.32803998 -0.35391774  0.61441445  0.57356854  0.24626412]
[ 0.45681372  0.17336728 -0.09373732 -0.04755792 -0.01471638]]




In [4]:

load = mat('0;0;0;2;-1')
qload = evecs.T*load ; Gamma = diag(asarray(qload).flatten())
print Gamma




[[ 0.19926624  0.          0.          0.          0.        ]
[ 0.         -0.88120277  0.          0.          0.        ]
[ 0.          0.          1.32256622  0.          0.        ]
[ 0.          0.          0.          1.194695    0.        ]
[ 0.          0.          0.          0.          0.50724461]]




In [5]:

R = mass*evecs*Gamma
print R




[[ 0.0076942   0.23355764 -0.69597406  0.72473136 -0.27000914]
[ 0.02204163  0.48367129 -0.67578503 -0.15422707  0.32429917]
[ 0.04177387  0.54535261  0.29379099 -0.6291333  -0.25178417]
[ 0.06536729  0.3118733   0.81260379  0.68523947  0.12491615]
[ 0.36411021 -0.61108692 -0.49589523 -0.22726885 -0.02985922]]




In [6]:

print [(R[i,:]).sum() for i in range(5)]
print [(R[:,i]).sum() for i in range(5)]




[-1.6653345369377348e-16, -1.1102230246251565e-16, -1.1102230246251565e-16, 1.9999999999999996, -0.99999999999999944]
[0.50098720202849911, 0.96336792762634316, -0.76125953450619877, 0.39934161163840598, -0.10243720678704951]




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...'




Now, the modal flexibility matrices, F_i = \psi*\psi^T/\omega^2,

F_1 =
[[ 0.067348  0.192931  0.365648  0.572162  0.796767]
[ 0.192931  0.552691  1.047475  1.639077  2.282505]
[ 0.365648  1.047475  1.985201  3.10642   4.325861]
[ 0.572162  1.639077  3.10642   4.860893  6.769061]
[ 0.796767  2.282505  4.325861  6.769061  9.426289]]
F_2 =
[[ 0.047616  0.098607  0.111182  0.063582 -0.031146]
[ 0.098607  0.204204  0.230245  0.131671 -0.0645  ]
[ 0.111182  0.230245  0.259608  0.148463 -0.072725]
[ 0.063582  0.131671  0.148463  0.084902 -0.04159 ]
[-0.031146 -0.0645   -0.072725 -0.04159   0.020373]]
F_3 =
[[ 0.024193  0.023492 -0.010213 -0.028248  0.00431 ]
[ 0.023492  0.02281  -0.009916 -0.027428  0.004185]
[-0.010213 -0.009916  0.004311  0.011924 -0.001819]
[-0.028248 -0.027428  0.011924  0.032981 -0.005032]
[ 0.00431   0.004185 -0.001819 -0.005032  0.000768]]
F_4 =
[[ 0.010428 -0.002219 -0.009052  0.00986  -0.000818]
[-0.002219  0.000472  0.001926 -0.002098  0.000174]
[-0.009052  0.001926  0.007858 -0.008559  0.00071 ]
[ 0.00986  -0.002098 -0.008559  0.009322 -0.000773]
[-0.000818  0.000174  0.00071  -0.000773  0.000064]]
F_5 =
[[ 0.004261 -0.005118  0.003974 -0.001971  0.000118]
[-0.005118  0.006147 -0.004773  0.002368 -0.000141]
[ 0.003974 -0.004773  0.003706 -0.001838  0.00011 ]
[-0.001971  0.002368 -0.001838  0.000912 -0.000055]
[ 0.000118 -0.000141  0.00011  -0.000055  0.000003]]

The following matrix, F-sum(F_i), should be close to zero...




In [9]:

print f
print "... and so it seems."




[[-0. -0. -0. -0. -0.]
[-0. -0. -0. -0. -0.]
[-0. -0. -0. -0. -0.]
[-0. -0. -0. -0. -0.]
[-0. -0. -0. -0. -0.]]
... 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"




The correction matrix for modes 3,4,5
[[ 0.038883  0.016154 -0.015291 -0.020359  0.00361 ]
[ 0.016154  0.02943  -0.012763 -0.027159  0.004217]
[-0.015291 -0.012763  0.015875  0.001527 -0.001   ]
[-0.020359 -0.027159  0.001527  0.043216 -0.005859]
[ 0.00361   0.004217 -0.001    -0.005859  0.000835]]

x = \psi_1 q_1 + \psi_2 q_2 + ((F-F_1-F_2)*r) f(t)
^^^^^^^^^^^^^^^
constant



let's try $p(t) = \boldsymbol{r} \\,\sin(t)$, that is $\omega=1.0$



In [ ]:




In [ ]: