```
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())
```




<style>
@font-face {
    font-family: 'Charis SIL';
    src: url('fonts/CharisSILEur-R.woff') format('woff');
}
@font-face {
    font-family: 'Architects Daughter';
    font-style: normal;
    src: url(https://fonts.gstatic.com/s/architectsdaughter/v6/RXTgOOQ9AAtaVOHxx0IUBM3t7GjCYufj5TXV5VnA2p8.woff2) format('woff2');
    unicode-range: U+0000-00FF, U+0131, U+0152-0153, U+02C6, U+02DA, U+02DC, U+2000-206F, U+2074, U+20AC, U+2212, U+2215, U+E0FF, U+EFFD, U+F000;
}
.prompt{display: None;}
div.cell{margin: auto;width:900px;}

div #notebook { /* centre the content */
    background: #fff9f0;
    margin: auto;
    padding: 0em;
    padding-top: 1em;
}
div #notebook_container {width: 960px!important;}
#notebook li { /* More space between bullet points */
    margin-top:0.2em;
}

div.text_cell_render{
    font-family: "Charis SIL", Cambria,  serif;
    line-height: 155%;
    font-size: 130%;
}
.CodeMirror{
    font-family: Consolas, monospace;
    font-size: 140%;
    background-color:#fffcf6!important;
}

.output_area {
    font-family: Consolas,monospace;font-size: 120%;
    background-color:#fffcf6!important;
    margin-top:0.8em;
}

div.output_latex{overflow:hidden}

h1 {
    font-family: 'Architects Daughter', serif;
    font-size:   48pt!important;
    text-align:  center;
    text-shadow: 4px 4px 4px #aaa;
    padding-bottom: 48pt;
}

.warning{color: rgb( 240, 20, 20 )}
</style>





# Rayleigh Quotient

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


```
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,


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

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


```
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()
```

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


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

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,


```
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()
```

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


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

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


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

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


```
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
```

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


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

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


```
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]]


921611d1f66d85dad852835c25a569d9f5d9cf8a
ipynb
Jupyter Notebook
dati_2015/ha04/03_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2015/ha04/03_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
dati_2015/ha04/03_Rayleigh.ipynb
shishitao/boffi_dynamics
365f16d047fb2dbfc21a2874790f8bef563e0947
MIT
2019-06-23T12:32:39.000Z
2021-08-15T18:33:55.000Z
Qwen/Qwen-72B
1. YES
2. YES