Rayleigh Quotient

Preliminaries

Two helper functions

def cumsum(lst):
    res = [0]
    for n in lst: 
        res.append(res[-1]+n)
    return res[1:]
def diff(lst):
    return [a-b for a, b in zip(lst, [0]+lst[:-1])]

Data

The storey stiffnesses, the storey masses and the assumed initial shape of the free vibrations.

ks = [5, 4, 2]
ms = [3, 3, 2]
xs = [1, 1, 1]

Direct Approach

Or, we use the storey stiffnesses to compute interstorey deflections, etc etc

Max Energies, R00

We compute twice the maximum strain and kinetic energies, omitting all the scaling factors from the expressions, but we use (formally) the scaling factors when we print the values. Eventually we compute the 1st approximation to the 1st eigenvalue.

V0 = sum(k*dx*dx for k, dx in zip(ks, diff(xs)))
T0 = sum(m*x*x for m, x in zip(ms, xs))
R00 = V0/T0

print('2V0 =', '%.3f'%V0, 'Z^2 k') 
print('2T0 =', '%.3f'%T0, 'Z^2 m w^2')
print('R00 =', '%.3f'%R00, 'k/m')

2V0 = 5.000 Z^2 k
2T0 = 8.000 Z^2 m w^2
R00 = 0.625 k/m

Better Strain Energy, R01

The inertial forces, the story shears using cumsum, the relative storey deflections and eventually, using again cumsum, the displacements due to the inertial forces.

fi = [m*x for m, x in zip(ms, xs)]
print('Inertial forces:         ', fi, 'Z m w^2')
ss = cumsum(fi[::-1])[::-1]
print('Storey shears:           ', ss, 'Z m w^2')
dx = [sh/k for sh, k in zip(ss, ks)]
print('Inter-storey deflections:', dx, 'Z m w^2 / k')
x1 = cumsum(dx)
print('Displacements:           ', x1, 'Z m w^2 / k')

Inertial forces:          [3, 3, 2] Z m w^2
Storey shears:            [8, 5, 2] Z m w^2
Inter-storey deflections: [1.6, 1.25, 1.0] Z m w^2 / k
Displacements:            [1.6, 2.85, 3.85] Z m w^2 / k

The strain energy and the new approximation to $\omega^2$:

V1 = sum(f*x for f, x in zip(fi, x1) )
print('2V1 =', '%.3f'%V1, 'Z^2 m^2 w^4 / k')
print('R01 =', '%.3f'%(T0/V1), 'k/m')

2V1 = 21.050 Z^2 m^2 w^4 / k
R01 = 0.380 k/m

Better Kinetic Energy, R11

T1 = sum(m*x*x for m, x in zip(ms, x1))
print('2T1 =', '%.3f'%T1, 'Z^2 m^3 w^6 / k^2')
print('R11 =', '%.3f'%(V1/T1), 'k/m')

2T1 = 61.693 Z^2 m^3 w^6 / k^2
R11 = 0.341 k/m

Matrix Approach

We define the structural matrices, the trial vector and compute the max energies in terms of double matrix products:

from numpy import array
M  = array(((3, 0, 0),
            (0, 3, 0),
            (0, 0, 2)), dtype='float')
K  = array(((9, -4, 0), (-4, 6, -2), (0, -2, 2)), dtype='float')
F  = array(((4, 4, 4), (4, 9, 9), (4, 9, 19)), dtype='float')/20
x0 = array((1,1,1), dtype='float')
V0 = x0@K@x0
T0 = x0@M@x0
print('2V0 =', '%.3f'%V0, 'Z^2 k') 
print('2T0 =', '%.3f'%T0, 'Z^2 m w^2')
R00 = V0/T0

print('R00 =', '%.3f'%R00, 'k/m')

2V0 = 5.000 Z^2 k
2T0 = 8.000 Z^2 m w^2
R00 = 0.625 k/m

It's much easier to compute the displacements due to the inertial forces...

fi = M@x0
x1 = F@fi
V1 = fi@x1
print()
print('2V1 =', '%.3f'%V1, 'Z^2 m^2 w^4 / k') 
print('R01 =', '%.3f'%(T0/V1), 'k/m')

2V1 = 21.050 Z^2 m^2 w^4 / k
R01 = 0.380 k/m

and the velocities and the kinetic energy...

T1 = x1@M@x1
print()
print('2T1 =', '%.3f'%T1, 'Z^2 m^3 w^6 / k^2')
print('R11 =', '%.3f'%(V1/T1), 'k/m')

2T1 = 61.693 Z^2 m^3 w^6 / k^2
R11 = 0.341 k/m

Cheat Check

Let's use a dedicated library function to check our best approximation.

from scipy.linalg import eigh
eigh(K, M)

(array([ 0.33333333,  1.66666667,  4.        ]),
 array([[-0.17407766,  0.3086067 , -0.45584231],
        [-0.34815531,  0.3086067 ,  0.34188173],
        [-0.52223297, -0.46291005, -0.11396058]]))

Not too bad indeed.