<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 2 DOF System
The data that describes our system and its loading
Step1:  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$.
Step2:  The @ operator stands, in this context, for matrix multiplication.
We print the computed values.
Step3:  Modal Response
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).$$
Step4:  The definition of time vector is a bit complicated...
Step5:  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.
<ASSISTANT_TASK:>
Python Code:
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')
Explanation: 2 DOF System
The data that describes our system and its loading:
End of explanation
evals, Psi = eigh(K, M)
Mstar = Psi.T@M@Psi
Kstar = Psi.T@K@Psi
pstar = Psi.T@p
Explanation: 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$.
End of explanation
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')
Explanation: The @ operator stands, in this context, for matrix multiplication.
We print the computed values.
End of explanation
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))
Explanation: Modal Response
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).$$
End of explanation
curves = plt.plot(t,q)
plt.legend(curves,['q1', 'q2'])
plt.title('Modal Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$q_i/\Delta_{st}$');
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}$');
Explanation: The definition of time vector is a bit complicated...
End of explanation
%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'$$'))
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()
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()
from IPython.display import HTML
HTML(open('00_custom.css').read())
Explanation: 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.
End of explanation
<END_TASK>