2 DOF System

The data that describes our system and its loading:



In [3]:

    
M = np.array(((2.0, 0.0), ( 0.0, 1.0)))
K = np.array(((3.0,-2.0), (-2.0, 2.0)))
p = np.array(( 0.0, 1.0)); w = 2.0

print_mat(M, pre='\\boldsymbol{M}=m\\,', fmt='%d')
print_mat(K, pre='\\boldsymbol{K}=k\\,', fmt='%d')
print_mat(p[:,None], pre=r'\boldsymbol{p}(t) = p_0\,', fmt='%d',
          post='\\sin(%d\\omega_0t)'%w, mt='B')









    





$$\boldsymbol{M}=m\,\begin{bmatrix}2&0\\0&1\end{bmatrix}$$






    





$$\boldsymbol{K}=k\,\begin{bmatrix}3&-2\\-2&2\end{bmatrix}$$






    





$$\boldsymbol{p}(t) = p_0\,\begin{Bmatrix}0\\1\end{Bmatrix}\sin(2\omega_0t)$$

Computing the eigenvalues and the eigenvectors

Using the eigh library function, specialized for Hermitian (symmetric) matrices.

We compute also the starred matrices. Note that the eigenvectors returned by eigh are normalized with respect to M, i.e., $M_i=m$ for every $i$.



In [4]:

    
evals, Psi = eigh(K, M)
Mstar = Psi.T@M@Psi
Kstar = Psi.T@K@Psi
pstar = Psi.T@p

The @ operator stands, in this context, for matrix multiplication.

We print the computed values.



In [5]:

    
print_mat(evals[None,:], mt='p', pre=r'\omega^2_i=\omega^2_o\,')
print_mat(Psi, pre=r'\boldsymbol{\Psi}=')
print_mat(Mstar, pre=r'\boldsymbol{M}^\star=m\,')
print_mat(Kstar, pre=r'\boldsymbol{K}^\star=k\,')
print_mat(pstar[:,None], pre=r'\boldsymbol{p}^\star=p_o\,', mt='B')









    





$$\omega^2_i=\omega^2_o\,\begin{pmatrix}0.313859&3.186141\end{pmatrix}$$






    





$$\boldsymbol{\Psi}=\begin{bmatrix}-0.541774&-0.454401\\-0.642621&0.766185\end{bmatrix}$$






    





$$\boldsymbol{M}^\star=m\,\begin{bmatrix}1.000000&-0.000000\\-0.000000&1.000000\end{bmatrix}$$






    





$$\boldsymbol{K}^\star=k\,\begin{bmatrix}0.313859&0.000000\\0.000000&3.186141\end{bmatrix}$$






    





$$\boldsymbol{p}^\star=p_o\,\begin{Bmatrix}-0.642621\\0.766185\end{Bmatrix}$$

The modal equation of motion, divided by $m$ $$ \ddot q_i + \Lambda^2_i\omega_0^2 q_i = \frac{p_0}{k}\frac{k}{m}\,p^\star_i \sin2\omega_0t.$$

The particular integral $\xi_i(t)$ is given by

$$\xi_i = \Delta_{st} \frac{\omega_0^2}{(\Lambda^2_i-4)\omega_0^2}p^\star_i\sin2\omega_0t.$$

Eventually the modal response, computed in the hypotesis of initial rest conditions is

$$q_i = \frac{\Delta_{st}}{\Lambda^2_i-4}p^\star_i(\sin2\omega_0t-\beta\sin\Lambda_i\omega_0t).$$



In [6]:

    
L  = np.sqrt(evals)
DAF = 1.0/(L**2-w**2)
beta = w/L
t = np.linspace(0,60,601)[:,None]
q = pstar*DAF*(np.sin(w*t)-beta*np.sin(L*t))

The definition of time vector is a bit complicated...



In [7]:

    
curves = plt.plot(t,q)
plt.legend(curves,['q1', 'q2'])
plt.title('Modal Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$q_i/\Delta_{st}$');



In [8]:

    
x = (Psi@q.T).T

curves = plt.plot(t, x)
plt.legend(curves,['x1', 'x2'])
plt.title('Structural Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$X_i/\Delta_{st}$');

The following code cell (that is executed before any other code cell by the notebook) loads libraries (or functions from libraries) and determines the style of plots and of the notebook itself. Besides the cell defines a function to format conveniently matrices and vectors.



In [2]:

    
%matplotlib inline

import matplotlib.pyplot as plt
plt.style.use(['fivethirtyeight', '00_mplrc'])

import numpy as np ; from scipy.linalg import eigh
np.set_printoptions(suppress=False, linewidth=120)

from IPython.display import Latex

def print_mat(mat, pre='', post='', fmt='%.6f', mt='b'):
    display(Latex(
        r'$$' + pre + r'\begin{%smatrix}'%mt +
        r'\\'.join('&'.join(fmt%x for x in row) for row in mat) +
        r'\end{%smatrix}'%mt + post + r'$$'))



In [9]:

    
import sympy as sy
sy.init_printing(use_latex=1)
o = sy.symbols('Omega')
display(o)
sM = sy.Matrix(((2,0,),(0,1)))
sK = sy.Matrix(((3, -2),(-2,2)))
KooM = sK - o*sM
iKooM = KooM.inv()
sp = sy.Matrix(((0,),(1,)))
a,b=(iKooM*sp)
a.expand().simplify(), b.expand().simplify()









    





$$\Omega$$






    Out[9]:





$$\left ( \frac{2}{2 \Omega^{2} - 7 \Omega + 2}, \quad \frac{- 2 \Omega + 3}{2 \Omega^{2} - 7 \Omega + 2}\right )$$



In [10]:

    
with plt.style.context('classic'):
    plot = sy.plot(a, b, (o, 0, 5), ylim=(-10,10), 
                   show=False, adaptive=False, nb_of_points=601)
    plot[0].line_color = 'black'; plot[0].label = '$x_1$'
    plot[1].line_color = 'red'  ; plot[1].label = '$x_2$'
    plot.xlabel = '$\\beta^2$'; plot.ylabel = r'$x_{i,\mathrm{ss}}/\Delta_\mathrm{st}$'
    plot.legend = True
    plot.show()



In [11]:

    
from IPython.display import HTML
HTML(open('00_custom.css').read())









    Out[11]:



In [ ]: