<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <a id='top'></a>
Complex vibration modes
Complex vibration modes arise in experimental research and numerical simulations with non proportional damping. In such cases the eigenproblem is nonlinear, thus requiring additional work in order to obtain the vibration frequencies and mode shapes of the dynamic system.
The solution to the free vibration response of the dynamic equilibrium equation involves some form of linearization of the quadratic eigenvalue problem it yields. For more information on this subject see for example this report.
Table of contents
Preamble
Dynamic equilibrium equation
State space formulation
Dynamic system setup
Undamped system
Proportionally damped system
Non proportionally damped system
Conclusions
Odds and ends
Preamble
We will start by setting up the computational environment for this notebook. Furthermore, we will need numpy and scipy for the numerical simulations and matplotlib for the plots
Step1:  We will also need a couple of specific modules and a litle "IPython magic" to show the plots
Step2:  Back to top
Dynamic equilibrium equation
In structural dynamics the second order differential dynamic equilibrium equation can be written in terms of generalized coordinates (d[isplacement]) and their first (v[elocity]) and second (a[cceleration]) time derivatives
Step3:  We will also setup two damping matrices, one proportional to the mass and stiffness matrices (C1) and the other non proportional (C2)
Step4:  Back to top
Undamped system
In the undamped system the damping matrix is all zeros and therefore the eigenproblem is a linear one as it involves only the mass and stiffness matrices
Step5:  The angular frequencies are computed as the square root of the eigenvalues
Step6:  The modal vectors, the columns of the modal matrix, have unit norm
Step7:  Contrary to what is normaly done, we will visualize the modal vectors in a polar plot of the corresponding amplitudes and angles of equivalent complex values
Step8:  Back to top
Proportionally damped system
This damping matrix is orthogonal because the mass and stiffness matrices are also orthogonal
Step9:  The system and input matrices are the following
Step10:  The eigenanalysis yields the eigenvalues and eigenvectors
Step11:  As we can see, the eigenvalues come in complex conjugate pairs. Let us take only the ones in the upper half-plane
Step12:  These complex eigenvalues can be decomposed into angular frequency and damping coefficient
Step13:  The columns of the modal matrix, the modal vectors, also come in conjugate pairs, each vector having unit norm
Step14:  Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices
Step15:  We will visualize again the complex valued modal vectors with a polar plot of the corresponding amplitudes and angles
Step16:  Back to top
Non proportionally damped system
In non proportionally damped systems the damping matrix is not proportional neither to the mass matrix nor the stiffness matrix.
Non proportinal damping carries the fact that the damping matrix is not orthogonal anymore
Step17:  The system and input matrices are the following
Step18:  The eigenanalysis yields the eigenvalues and eigenvectors of the system matrix
Step19:  As we can see, the eigenvalues come in complex conjugate pairs. Again, let us take only the ones in the upper half-plane
Step20:  These complex eigenvalues can be decomposed into angular frequency and damping coefficient much like in the propotional damping case
Step21:  Again, the columns of the modal matrix, the modal vectors, come in conjugate pairs, and each vector has unit norm
Step22:  Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices
Step23:  Once more we will visualize the complex valued modal vectors through a polar plot of the corresponding amplitudes and angles
<ASSISTANT_TASK:>
Python Code:
import sys
import numpy as np
import scipy as sp
import matplotlib as mpl
print('System: {}'.format(sys.version))
print('numpy version: {}'.format(np.__version__))
print('scipy version: {}'.format(sp.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
Explanation: <a id='top'></a>
Complex vibration modes
Complex vibration modes arise in experimental research and numerical simulations with non proportional damping. In such cases the eigenproblem is nonlinear, thus requiring additional work in order to obtain the vibration frequencies and mode shapes of the dynamic system.
The solution to the free vibration response of the dynamic equilibrium equation involves some form of linearization of the quadratic eigenvalue problem it yields. For more information on this subject see for example this report.
Table of contents
Preamble
Dynamic equilibrium equation
State space formulation
Dynamic system setup
Undamped system
Proportionally damped system
Non proportionally damped system
Conclusions
Odds and ends
Preamble
We will start by setting up the computational environment for this notebook. Furthermore, we will need numpy and scipy for the numerical simulations and matplotlib for the plots:
End of explanation
from numpy import linalg as LA
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: We will also need a couple of specific modules and a litle "IPython magic" to show the plots:
End of explanation
MM = np.matrix(np.diag([1., 2.]))
print(MM)
KK = np.matrix([[20., -10.], [-10., 10.]])
print(KK)
Explanation: Back to top
Dynamic equilibrium equation
In structural dynamics the second order differential dynamic equilibrium equation can be written in terms of generalized coordinates (d[isplacement]) and their first (v[elocity]) and second (a[cceleration]) time derivatives:
\begin{equation}
\mathbf{M} \times \mathbf{a(t)} + \mathbf{C} \times \mathbf{v(t)} + \mathbf{K} \times \mathbf{d(t)} = \mathbf{F(t)}
\end{equation}
where:
$\mathbf{M}$ is the mass matrix
$\mathbf{C}$ is the damping matrix
$\mathbf{K}$ is the stiffness matrix
$\mathbf{a(t)}$ is the acceleration vector
$\mathbf{v(t)}$ is the velocity vector
$\mathbf{d(t)}$ is the displacement vector
$\mathbf{F(t)}$ is the force input vector  
Considering a dynamic system with $NDOF$ is the number of generalized degrees of freedom, the vectors will have dimensions of $NDOF \times 1$ and the matrices $NDOF \times NDOF$.
When the system is undamped, the damping matrix will be null. In this case the eigenproblem is linear:
\begin{equation}
\left[ -\mathbf{M} \times \mathbf{\omega^2} + \mathbf{K} \right] \times \mathbf{v} = \mathbf{0}
\end{equation}
In a proportionally damped system, the damping matrix is proportional to the mass and stiffness matrices:
\begin{equation}
\mathbf{C} = \alpha \times \mathbf{M} + \beta \times \mathbf{K}
\end{equation}
where $\alpha$ and $\beta$ are mass and stiffness proportionality coefficients. These are typically very small positive numbers. The resulting eigenproblem is still linear because the damping matrix can be decomposed by the modal vectors.
In a system with non proportional damping, the damping matrix will not be proportional neither to the mass nor to the stiffness matrices.
Back to top
Dynamic system setup
In this example we will use the folowing mass and stiffness matrices:
End of explanation
C1 = 0.1*MM+0.04*KK
print(C1)
C2 = np.matrix([[0.1, 0.2], [0.2, 0.2]])
print(C2)
Explanation: We will also setup two damping matrices, one proportional to the mass and stiffness matrices (C1) and the other non proportional (C2):
End of explanation
W2, F1 = LA.eig(LA.solve(MM,KK)) # eigenanalysis
ix = np.argsort(np.absolute(W2)) # sort eigenvalues in ascending order
W2 = W2[ix] # sorted eigenvalues
F1 = F1[:,ix] # sorted eigenvectors
print(np.round_(W2, 4))
print(np.round_(F1, 4))
Explanation: Back to top
Undamped system
In the undamped system the damping matrix is all zeros and therefore the eigenproblem is a linear one as it involves only the mass and stiffness matrices:
End of explanation
print(np.sqrt(W2))
Explanation: The angular frequencies are computed as the square root of the eigenvalues:
End of explanation
print(LA.norm(F1, axis=0))
Explanation: The modal vectors, the columns of the modal matrix, have unit norm:
End of explanation
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
    ax[mode].set_title('Mode #{}'.format(mode+1))
    for dof in range(2):
        r = np.array([0, np.absolute(F1[dof,mode])])
        t = np.array([0, np.angle(F1[dof,mode])])
        ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
Explanation: Contrary to what is normaly done, we will visualize the modal vectors in a polar plot of the corresponding amplitudes and angles of equivalent complex values:
End of explanation
print(np.round_(F1.T*C1*F1, 4))
Explanation: Back to top
Proportionally damped system
This damping matrix is orthogonal because the mass and stiffness matrices are also orthogonal:
End of explanation
A = np.bmat([[np.zeros_like(MM), MM], [MM, C1]])
print(A)
B = np.bmat([[MM, np.zeros_like(MM)], [np.zeros_like(MM), -KK]])
print(B)
Explanation: The system and input matrices are the following:
End of explanation
w1, v1 = LA.eig(LA.solve(A,B))
ix = np.argsort(np.absolute(w1))
w1 = w1[ix]
v1 = v1[:,ix]
print(np.round_(w1, 4))
print(np.round_(v1, 4))
Explanation: The eigenanalysis yields the eigenvalues and eigenvectors:
End of explanation
print(np.round_(w1[::2], 4))
Explanation: As we can see, the eigenvalues come in complex conjugate pairs. Let us take only the ones in the upper half-plane:
End of explanation
zw = -w1.real # damping coefficient time angular frequency
wD = w1.imag # damped angular frequency
zn = 1./np.sqrt(1.+(wD/-zw)**2) # the minus sign is formally correct!
wn = zw/zn # undamped angular frequency
print('Angular frequency: {}'.format(wn[[0,2]]))
print('Damping coefficient: {}'.format(zn[[0,2]]))
Explanation: These complex eigenvalues can be decomposed into angular frequency and damping coefficient:
End of explanation
print(LA.norm(v1[:,::2], axis=0))
Explanation: The columns of the modal matrix, the modal vectors, also come in conjugate pairs, each vector having unit norm:
End of explanation
AA = v1[:2,[0,2]]
AB = AA.conjugate()
BA = np.multiply(AA,w1[[0,2]])
BB = BA.conjugate()
v1_new = np.bmat([[AA, AB], [BA, BB]])
print(np.round_(v1_new[:,[0,2,1,3]], 4))
Explanation: Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices:
End of explanation
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
    ax[mode].set_title('Mode #{}'.format(mode+1))
    for dof in range(2):
        r = np.array([0, np.absolute(v1[dof,2*mode])])
        t = np.array([0, np.angle(v1[dof,2*mode])])
        ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
Explanation: We will visualize again the complex valued modal vectors with a polar plot of the corresponding amplitudes and angles:
End of explanation
print(np.round_(F1.T*C2*F1, 4))
Explanation: Back to top
Non proportionally damped system
In non proportionally damped systems the damping matrix is not proportional neither to the mass matrix nor the stiffness matrix.
Non proportinal damping carries the fact that the damping matrix is not orthogonal anymore:
End of explanation
A = np.bmat([[np.zeros_like(MM), MM], [MM, C2]])
print(A)
B = np.bmat([[MM, np.zeros_like(MM)], [np.zeros_like(MM), -KK]])
print(B)
Explanation: The system and input matrices are the following:
End of explanation
w2, v2 = LA.eig(LA.solve(A,B))
ix = np.argsort(np.absolute(w2))
w2 = w2[ix]
v2 = v2[:,ix]
print(np.round_(w2, 4))
print(np.round_(v2, 4))
Explanation: The eigenanalysis yields the eigenvalues and eigenvectors of the system matrix:
End of explanation
print(np.round_(w2[[0,2]], 4))
Explanation: As we can see, the eigenvalues come in complex conjugate pairs. Again, let us take only the ones in the upper half-plane:
End of explanation
zw = -w2.real # damping coefficient times angular frequency
wD = w2.imag # damped angular frequency
zn = 1./np.sqrt(1.+(wD/-zw)**2) # the minus sign is formally correct!
wn = zw/zn # undamped angular frequency
print('Angular frequency: {}'.format(wn[[0,2]]))
print('Damping coefficient: {}'.format(zn[[0,2]]))
Explanation: These complex eigenvalues can be decomposed into angular frequency and damping coefficient much like in the propotional damping case:
End of explanation
print(LA.norm(v2[:,[0,2]], axis=0))
Explanation: Again, the columns of the modal matrix, the modal vectors, come in conjugate pairs, and each vector has unit norm:
End of explanation
AA = v2[:2,[0,2]]
AB = AA.conjugate()
BA = np.multiply(AA,w2[[0,2]])
BB = BA.conjugate()
v2_new = np.bmat([[AA, AB], [BA, BB]])
print(np.round_(v2_new[:,[0,2,1,3]], 4))
Explanation: Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices:
End of explanation
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
    ax[mode].set_title('Mode #{}'.format(mode+1))
    for dof in range(2):
        r = np.array([0, np.absolute(v2[dof,2*mode])])
        t = np.array([0, np.angle(v2[dof,2*mode])])
        ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
Explanation: Once more we will visualize the complex valued modal vectors through a polar plot of the corresponding amplitudes and angles:
End of explanation
<END_TASK>