First, initialize our environment
In [1]:
### Plotting stuff
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.autolayout'] = True
plt.style.use(['fivethirtyeight', './00_mplrc'])
figd = '../figures/05_'
### math stuff
import numpy as np
np.set_printoptions(linewidth=120)
from scipy.linalg import eigh
### format the notebook, ancillary stuff
from IPython.display import display, Latex, HTML
display(HTML(open("00_custom.css", "r").read()))
Next, some handy functions have to be defined...
In [2]:
### A bunch of functions to help with repetitive tasks
def pmat(mat, fmt="%f+8.3", delim="b"):
"pmat prints an array to be 'copy&paste'd into a LaTeX document"
print
print r'\begin{'+delim+'matrix}'
print '\\\\\n'.join(' & '.join(fmt%n for n in row) for row in mat)
print r'\end{'+delim+'matrix}'
print
##########
def xdot(*mats):
"xdot(m1,m2, ..., mn) computes a dot product of all the arguments"
return reduce(np.dot, mats)
##########
def norm(v):
"norm(v) normalizes a vector so that its norm is 1"
return v/np.sqrt(np.dot(v, v))
##########
def normalize(thearray, M):
"normalize(v, M) normalizes a vector so that the norm v.T*M*v is 1"
return thearray/np.sqrt(np.diag(thearray.T.dot(M.dot(thearray))))
### I've a couple of aux functions to have customized plots of eigenvectors
### in an external module
from eigenplot import plot_evec, plot_ritz
Starting from a flist with the number of each floor, we generate the masses of each floor and the stiffnesses below of each floor, and add an additional, fictitiuos zero stiffness to our list to represent the stiffness of the non-existing post-ultimate floor...
From these lists, using np.diag it's easy to construct the structural matrices.
Eventually we print our structural matrices, using a normalization factor to keep the width of the printing tolerable.
In [3]:
n_floor = range(1,10) # --> (1, ..., 9)
m_floor = [405000. - floor_n*5000. for floor_n in n_floor]
k_floor = [(840-10*(floor_n-1)**2)*1E6 for floor_n in n_floor]
k_floor.append(0.0)
M = np.diag(m_floor)
K = (np.diag(k_floor[:-1]) +
np.diag(k_floor[+1:]) -
np.diag(k_floor[1:-1],k=+1) -
np.diag(k_floor[1:-1],k=-1))
print ";".join("%d"%(m/1000) for m in m_floor)
print ";".join("%d"%(m/1E6) for m in k_floor)
print 'mass matrix M/1000\n', M/1000
print 'stiffness matrix K/1E6\n', K/1E6
pmat(M/1000, fmt='%d')
pmat(K/1E6, fmt='%d')
Instead of computing the exact eigenvalues at the end, we do it immediately using scipy.linalg.eigh (here used as eigh), so that we can do our comparisons while we compute new approximations.
In [4]:
evls, evcs = eigh(K,M,eigvals=(0,2))
Do we want to plot the eigenvectors? Definitely yes!
In [5]:
plot_evec(evcs, 'The "real" lowest 3 eigenvectors',
save=figd+'real_ev.pdf')
print "The eigenvalues, the natural periods"
print evls
print np.array([2*np.pi/np.sqrt(evl) for evl in evls])
print
print "The eigenvectors"
print evcs
In [6]:
phi0 = np.array((( 0.33, 0.67, 1.00, 1.33, 1.67, 2.00, 2.33, 2.67, 3.00),
(-1.00,-2.00,-3.00,-3.00,-2.00,-1.00, 0.00, 1.00, 3.00),
( 1.00, 2.00, 2.00, 0.00,-2.00,-2.00, 0.00, 1.00, 3.00))).T
plot_evec(phi0, '$\\Phi_0$, not normalized', save=figd+'phi0_nnorm.pdf')
pmat(phi0.T, fmt='%+05.2f')
phi0 = normalize(phi0, M)
plot_evec(phi0, '$\\Phi_0$, normalized', save=figd+'phi0_norm.pdf')
Mz = np.dot(phi0.T,np.dot(M,phi0))
Kz = np.dot(phi0.T,np.dot(K,phi0))
zvl, zvc = eigh(Kz, Mz)
zvc*= np.sign(phi0.dot(zvc)[-1])
print evls
print zvl
print zvc
plot_ritz(zvc, save=figd+'zvc0.pdf')
plot_evec(np.dot(phi0,zvc), 'The lowest 3 eigenvectors, 1st estimate',
save=figd+'psi_est0.pdf')
pmat(zvc, fmt='%+9.6f')
First, the new base $\Phi_1 = K^{-1}M\Phi_0z$
In [7]:
phi1 = xdot(np.linalg.inv(K),M,phi0,zvc)
phi1 = normalize(phi1, M)
display(HTML("<h3>The new base</h3>"))
plot_evec(phi1,'The new Ritz base $\\Phi_1$, normalized',
save=figd+'phi1_norm.pdf')
In [8]:
Mz = np.dot(phi1.T,np.dot(M,phi1))
Kz = np.dot(phi1.T,np.dot(K,phi1))
zvl, zvc = eigh(Kz, Mz)
zvc*= np.sign(phi1.dot(zvc)[-1])
print evls
print zvl
print zvc
pmat(zvc, fmt='%+9.6f')
plot_ritz(zvc, save=figd+'zvc1.pdf')
plot_evec(np.dot(phi1,zvc), 'The lowest 3 eigenvectors, 2nd estimate',
save=figd+'psi_est1.pdf')
In [9]:
plot_evec(np.dot(phi1,zvc)-evcs,
title='Differences, "real" - estimate', save=figd+'diff.pdf')
In [9]: