Some imports
In [1]:
from scipy import * ; from scipy.linalg import eigh
One utility function, to print a matrix in a format recognized by typesetting program
In [2]:
def lp(data,name,fmt="%+10.4f",title=""):
delim={"mat":"b", "vet":"B", "det":"V", "norm":"v"}
if title: print "% ----- "+title+" -----"
print "\\begin{"+delim[name]+"matrix}"
print "\\\\\n".join(["&".join(map(lambda x: fmt%(x,),line))
for line in asarray(data)])
print "\\end{"+delim[name]+"matrix}"
Here a simple function to compute the DRV's as a function of the NDOFs "n" and of the load vector "r"
In [3]:
def derritz(n,r):
#global T, phi
phi = mat(zeros((n,n)))
# y are the initial displacements
y = F*r
# normalization factor, note that M = I
b = sqrt((y.T*y)[0,0])
# normalize
y = y/b
# save the first DRV
phi[:,0] = y[:,0]
for i in range(1,n):
# remember that M = I
y = F*phi[:,i-1]
alpha = y.T*phi
for k in range(i):
y=y-alpha[0,k]*phi[:,k]
b=sqrt((y.T*y)[0,0])
phi[:,i]=y[:,0]/b
# remember that M = I
T=phi.T*F*phi
# there are very small terms outside the tridiagonal part
# this reduces T to a true tridiagonal matrix...
T=(diagflat(diag(T))+
diagflat(diag(T,1),1)+
diagflat(diag(T,-1),-1))
return phi, T
The problem data:
In [4]:
n=5
M=mat(eye(n))
K=mat(eye(n))*2
K[n-1,n-1]=1.
for i in range(n-1):
K[i+1,i]=-1.
K[i,i+1]=-1.
F=K.I
lp(M,'mat',title="Mass matrix (=I) M")
lp(K,'mat',title="Stiffness matrix K")
lp(F,'mat',title="Flexibility matrix F")
# evecs are normalized with respect to M
evals, evecs = eigh(K,M)
L=mat(diagflat(evals))
lp(L,'mat',title="Eigenvalues Matrix L")
lp(evecs,'mat',title="Eigenvectors matrix, \Psi")
Let's set a particular load vector, compute the DRVs and show the results
In [5]:
r = mat((0,0,0,-2.,1.)).T
lp(r.T,'vet',title="Load vector transposed r^T")
phi, T = derritz(n,r)
lp(phi,'mat',title="Derived Ritz Vectors matrix \Phi")
lp(T,'mat',title="Tridiagonal DRV matrix, T")
Compute the modal partecipation factors, for the eigenvectors and for the DRVs, and the decomposition in modal forces and DRV forces
In [6]:
Gamma=evecs.T*r ; Gammh=phi.T*r
f_m=M*evecs*diagflat(Gamma)
f_r=M*phi*diagflat(Gammh)
lp(Gamma.T,'vet',title="Modal partecipation factors")
lp(Gammh.T,'vet',title="DRV's partecipation factors")
lp(f_m,'mat',title="Modal forces matrix")
lp(f_r,'mat',title="DRV's forces matrix")
In [7]:
den=dot(r.T,r)[0,0] ; e_m = r ; e_r = r
for i in range(n):
e_m = e_m-f_m[:,i]
e_r = e_r-f_r[:,i]
print "%3d %10.7f %10.7f" % (
i+1,
dot(r.T,e_m)[0,0]/den,
dot(r.T,e_r)[0,0]/den )
In [7]: