

In [1]:

    
from sympy import *
init_printing(use_latex=True)
from IPython.display import display
from numpy import array
from scipy.linalg import eigh
from IPython.display import HTML
HTML(open("00_custom.css", "r").read())









    Out[1]:

Rayleigh Quotient

We are going to use a little symbolic algebra, hence we'll need the following base symbols:



In [2]:

    
x1, x2, x3, k, L, m = symbols('x1, x2, x3, k, L, m')

The single degrees of freedom have to be collected in a vector,



In [3]:

    
x = Matrix((x1,x2,x3))
x









    Out[3]:





$$\left[\begin{matrix}x_{1}\\x_{2}\\x_{3}\end{matrix}\right]$$

We want to compute the (double of the) strain energy, to do so we need the relative rotations, that we can compute if we know in advance the rotations of each individual bar...



In [9]:

    
phi1, phi2, phi3 = x1/L, (x2-x1)/L, (x3-x2)/L
phi21, phi32 = phi2-phi1, phi3-phi2
W2 = x2**2*k + phi21**2*k*L*L + phi32**2*k*L*L
print latex((W2/k).expand())
W2.expand()









    



5 x_{1}^{2} - 8 x_{1} x_{2} + 2 x_{1} x_{3} + 6 x_{2}^{2} - 4 x_{2} x_{3} + x_{3}^{2}






    Out[9]:





$$5 k x_{1}^{2} - 8 k x_{1} x_{2} + 2 k x_{1} x_{3} + 6 k x_{2}^{2} - 4 k x_{2} x_{3} + k x_{3}^{2}$$

By comparison with the previuous expression, we have the stiffness matrix



In [6]:

    
K = k*Matrix(((5,-4,1),(-4,6,-2),(1,-2,1)))
K









    Out[6]:





$$\left[\begin{matrix}5 k & - 4 k & k\\- 4 k & 6 k & - 2 k\\k & - 2 k & k\end{matrix}\right]$$

To compute the (double of the) kinetic energy, we have to compute the rotatory inertia and the displacements of the centers of the bars, as we have already computed the bars' rotations,



In [8]:

    
J = m*L*L/12
xg1, xg2, xg3 = x1/2, (x1+x2)/2, (x2+x3)/2
T2 = m*(xg1**2+xg2**2+xg3**2) + J*(phi1**2+phi2**2+phi3**2) 
print latex((T2/m).expand())
T2.expand()









    



\frac{2 x_{1}^{2}}{3} + \frac{x_{1} x_{2}}{3} + \frac{2 x_{2}^{2}}{3} + \frac{x_{2} x_{3}}{3} + \frac{x_{3}^{2}}{3}






    Out[8]:





$$\frac{2 m}{3} x_{1}^{2} + \frac{m x_{1}}{3} x_{2} + \frac{2 m}{3} x_{2}^{2} + \frac{m x_{2}}{3} x_{3} + \frac{m x_{3}^{2}}{3}$$

Again, by comparison with the kinetic energy we can derive the mass matrix,



In [11]:

    
M = m*Matrix(((4,1,0),(1,4,1),(0,1,2)))/6
M









    Out[11]:





$$\left[\begin{matrix}\frac{2 m}{3} & \frac{m}{6} & 0\\\frac{m}{6} & \frac{2 m}{3} & \frac{m}{6}\\0 & \frac{m}{6} & \frac{m}{3}\end{matrix}\right]$$

Having computed, as requested, the structural matrices we are ready for the Rayleigh estimates, but first we define a trial vector and an auxiliary matrix, that we'll use in the following to simplify our lengthy expressions



In [12]:

    
x0 = Matrix((1,2,4))
D = K.inv()*M

Here are the expression for the energies, normalized with respect to $Z_0^2/2$



In [13]:

    
w2 = 1
V0 = (x0.T*K*x0)[0] *     w2**0
T0 = (x0.T*M*x0)[0] *     w2**1
V1 = (x0.T*M*D*x0)[0] *   w2**2
T1 = (x0.T*M*D*D*x0)[0] * w2**3
V0, T0, V1, T1









    Out[13]:





$$\left ( 5 k, \quad 12 m, \quad \frac{355 m^{2}}{12 k}, \quad \frac{31523 m^{3}}{432 k^{2}}\right )$$

The Rayleigh estimates are as follows, fractions first and then floating point numbers



In [10]:

    
display(V0/T0, T0/V1, V1/T1)
map(N,(V0/T0, T0/V1, V1/T1))









    





$$\frac{5 k}{12 m}$$






    





$$\frac{144 k}{355 m}$$






    





$$\frac{12780 k}{31523 m}$$






    Out[10]:





$$\left [ \frac{0.416666666666667 k}{m}, \quad \frac{0.405633802816901 k}{m}, \quad \frac{0.40541826602798 k}{m}\right ]$$

Just to be sure, a check with the results of a library function...



In [11]:

    
for x in eigh(
    array(((5,-4,1),(-4,6,-2),(1,-2,1)))/1.,
    array(((4,1,0),(1,4,1),(0,1,2)))/6., eigvals=(0,0)): print x









    



[ 0.40540716]
[[ 0.2843138 ]
 [ 0.60585859]
 [ 1.11964382]]