In [1]:
from scipy import *
from IPython.display import display, HTML, Latex, SVG
from fractions import Fraction as frac
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
Out[1]:
In [2]:
display(SVG(filename='rigid.svg'))
For a single bar, $$\dot x(s) = \dot x_1 + \frac{\dot x_{i+1} - \dot x_i}{L}\,s$$ and the total kinetic energy is $$T=\frac12\int_0^L \frac{m}{L}\dot x^2(s)\,\text{d}s =\frac12\left( \frac13\dot x_i^2 + \frac13\dot x_i \dot x_{i+1} + \frac13\dot x_{i+i}^2 \right)\, m $$
We compute, bar by bar, the contribution to twice the kinetic energy. The velocities of the first and last bar have to be written in terms of the fictitious coordinates $x_0$ and $x_4$ where $x_0=x_4=0$
In [3]:
def index(i):
return i*2, i*2+1, i*2+2
# We use x0, x1, x2, x3, x4, with the understanding that x0=x4=0
# we compute twice the kinetic energy storing the coefficients
# of the quadratic form in a sequence of length 9 --- note
# that the cross terms are restricted to adjacent DOFs
# Initially T2 is zero
# x0x0, x0x1, x1x1, x1x2, x2x2, x2x3, x3x3, x3x4, x4x4
T2 = [ 0, 0, 0, 0, 0, 0, 0, 0, 0]
# then for each bar we add, in the correct positions, the contributions to T2
# due to the DOFs associated with the current bar
for i in (0,1,2,3):
i11,i12,i22=index(i)
T2[i11] += frac(1,3)
T2[i12] += frac(1,3)
T2[i22] += frac(1,3)
we print the coefficients, omitting the ones associated with either $x_0$ or $x_4$
In [4]:
print T2[2:-2]
With the help of a helper function, we display the quadratic form for $2T$ and also for $2V$, that's so easy to determine that I leave its determination to the reader
In [5]:
def Lf(fr):
n = fr.numerator
d = fr.denominator
if n>0:
return r' + \frac{%d}{%d} '%( n,d) if d>1 else '%d'%(n,)
elif n<0:
return r' - \frac{%d}{%d} '%(-n,d) if d>1 else '%d'%(n,)
else:
return ' 0 '
s = (r'$$2 T /m= ' +
Lf(T2[2])+'x_1^2' +
Lf(T2[3])+'x_1x_2'+
Lf(T2[4])+'x_2^2'+
Lf(T2[5])+'x_2x_3'
+Lf(T2[6])+'x_3^2$$')
display(Latex(s))
display(Latex(r'$$2V/k=\frac1{12}x_1^2+\frac1{12}x_1x_2+\frac1{12}x_2^2+x_3^2$$'))
Using again fraction, we construct the nass matrix and the stiffness matrix equating the double matrix products to our previous results.
In [6]:
six = frac(1,6)
M = matrix((
(4*six,six,0),
(six,4*six,six),
(0,six,4*six)
))
tw4 = frac(1,24)
one = frac(1)
K = matrix((
(2*tw4,tw4, frac()),
(tw4,2*tw4, frac()),
(frac(),frac(),one)))
Let's store the stiffness inverse (a bit of work to have integer coefficients).
In [7]:
F = matrix(map(int,ravel(K.I))).reshape((3,3))
Our initial guess for the shape vector
In [14]:
# u = matrix('-10;10;1')
u = matrix('1;-1;0')
Our initial guess for the Rayleigh estimate,
In [15]:
n, = ravel(u.T*K*u)
d, = ravel(u.T*M*u)
print n, d, '\t\tRoo =', n/d,'=',1.*n/d
for the second estimate, the previous denominator is the new numerator, and we have to compute a new denominator
In [16]:
n = d
d, = ravel(u.T*M*F*M*u)
print n, d, '\t\tRoi =', n/d,'=',1.*n/d
and then the final estimate that was requested
In [17]:
n = d
d, = ravel(u.T*M*F*M*F*M*u)
print n, d, '\t\tRii =', n/d,'=',1.*n/d
To see if ours result is good, we have to cheat...
In [12]:
from scipy.linalg import eigh
eval, evec = eigh(K,M,eigvals=(0,0))
print "1st eigenvalue = ", eval[0]
print "1st eigenvector =", ravel(evec)
The results are good, but I sort of cheated in the choice of the trial shape vector, as I already had an idea of the first eigenvector of the system.