
In the last lecture, we obtained the modal properties of the system by constructing the characteristic equation from $N$ simultaneous equations. Analysing the system by directly managing $N$ simultaneous equations in this way quickly becomes unwieldy for larger systems with many more degrees of freedom.

Therefore we cast our equations in matrix form and make use of the properties of matrices for larger systems. Our equations of motion in matrix form are,

$$
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\ddot{u}
\end{Bmatrix} +
\begin{bmatrix}
K
\end{bmatrix}
\begin{Bmatrix}
u
\end{Bmatrix}=0
\tag{1}
$$

where

$$
\begin{bmatrix}
M
\end{bmatrix} =
\begin{bmatrix}
m_1 & 0\\
0 & m_2
\end{bmatrix}
\tag{2}
$$

$$
\begin{bmatrix}
K
\end{bmatrix} =
\begin{bmatrix}
k_1+k_2 & -k_2\\
-k_2 & k_2
\end{bmatrix}
\tag{3}
$$

$$
\begin{Bmatrix}
u
\end{Bmatrix} =
\begin{Bmatrix}
u_1\\
u_2
\end{Bmatrix} \: \text{and}\:
\begin{Bmatrix}
\ddot{u}
\end{Bmatrix} =
\begin{Bmatrix}
\ddot{u}_1\\
\ddot{u}_2
\end{Bmatrix}
\tag{4}
$$

We again assume a harmonic solution where $u_1 = U_1\sin(\omega t)$ and $u_2 = U_2\sin(\omega t)$ or in matrix form,

$$
\begin{Bmatrix}
u
\end{Bmatrix} =
\begin{Bmatrix}
U
\end{Bmatrix} \sin(\omega t)
\tag{5}
$$

After differentiating we have,

$$
\tag{6} \begin{Bmatrix}
\ddot{u}
\end{Bmatrix} = -\omega^2
\begin{Bmatrix}
u
\end{Bmatrix}
$$

Substituting our assumed solutions back into the matrix equation of motion and factoring out $\sin(\omega t)$ yields,

$$
-\omega^2
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
U
\end{Bmatrix} +
\begin{bmatrix}
K
\end{bmatrix}
\begin{Bmatrix}
U
\end{Bmatrix} = 0
\tag{7}
$$

Rearranging slightly,

$$
\left(
\begin{bmatrix}
K
\end{bmatrix}-\omega^2
\begin{bmatrix}
M
\end{bmatrix}
\right)
\begin{Bmatrix}
U
\end{Bmatrix} = 0
\tag{8}
$$

Remembering that $\begin{Bmatrix}U\end{Bmatrix} = \begin{Bmatrix}U_1 U_2\end{Bmatrix}^T$; this equation is an example of a specific mathematical problem known as an **eigenvalue problem** or **eigenproblem**. By recognising this equation as an eigenvalue problem we know the following statement to be true,

$$
\left|
\begin{bmatrix}
K
\end{bmatrix}-\omega^2
\begin{bmatrix}
M
\end{bmatrix}
\right| = 0
\tag{9}
$$

i.e. the determinate of the matrix factor of $\{U\}$ must equal zero. Expanding this expression results in the characteristic polynomial equation we developed in the last lecture. Thus we have arrived at the same destination (the system's characteristic equation) by way of matrix notation.

### General Eigenvalue formulation

We can rewrite our matrix equations of motion to correspond with the classic representation or standard eigenproblem form (that you will likely see in a maths text). For completeness we'll demonstrate that here; the matrix equations of motion,

$$
\left(
\begin{bmatrix}
K
\end{bmatrix}-\omega^2
\begin{bmatrix}
M
\end{bmatrix}
\right)
\begin{Bmatrix}
U
\end{Bmatrix} = 0
\tag{10}
$$

can be rearranged to,

$$
\begin{bmatrix}
M
\end{bmatrix}^{-1}
\begin{bmatrix}
K
\end{bmatrix}
\begin{Bmatrix}
U
\end{Bmatrix} = \omega^2
\begin{Bmatrix}
U
\end{Bmatrix}
\tag{11}
$$

This corresponds to the standard eigenproblem,

$$
\begin{bmatrix}
A
\end{bmatrix}
\begin{Bmatrix}
U
\end{Bmatrix} = \lambda
\begin{Bmatrix}
U
\end{Bmatrix}
\tag{12}
$$

where $\begin{Bmatrix}U\end{Bmatrix}$ represents the eigenvector of matrix $\begin{bmatrix}A\end{bmatrix}$ and $\lambda$ represents the eigenvalue.

In words, we would say that an eigenvector is some vector that when premultiplied by a square matrix, yields another vector that is a scaled version of the original eigenvector. The scale factor is known as the eigenvalue. For a deeper investigation of the mathematics of eigenvalues and eigenvectors refer to your nearest maths textbook.

For our purposes in the study of structural dynamics, the key point is this, **the natural frequencies and corresponding mode shapes of the system can be obtained from the eigenvalues and eigenvectors respectively of the matrix $[M]^{-1}[K]$.**

So, once we have the mass and stiffness matrix of our system, we can construct the square matrix $[M]^{-1}[K]$ and find its eigenvalues and eigenvectors to yield the dynamic characteristics of the system. Calculation of eigenvalues and eigenvectors is a computationally expensive task for larger systems. Fortunately, others have done the hard work of developing methods and tools that we can deploy. So we won't be delving into the details of the calculation; we'll just reach for the tried and tested tools already in existence.

### Eigenvalues & Eigenvector solution using SciPy

The previous solution method relied on directly identifying the roots of the characteristic equation. As mentioned, this becomes more challenging for systems with many degrees of freedom. For larger systems, we can employ computational techniques to extract the eigenvalues and eigenvectors of the matrix $[M]^{-1}[K]$. We now know that these correspond to the system natural frequencies and mode shapes. We will make use of the **linalg** (linear algebra) collection of functions from the **SciPy** library to do the heavy lifting.

We'll start by confirming the natural frequencies and mode shapes for the 2DoF system analysed in the previous lecture. First we need to import the relevant tools from the SciPy library.,

```python
import scipy.linalg as la
```

Now we can construct the mass and stiffness matrices and pass the matrix $[M]^{-1}[K]$ to the `eig` function from the `linalg` toolbox.

```python
M = np.mat([[m1, 0],[0,m2]]) #Mass matrix
K = np.mat([[k1+k2, -k2],[-k2, k2]]) #Stiffness matrix
eValues, eVectors = la.eig(M.I*K) #Use library to extract eigenvalues and eigenvectors
```

We can now extract the modal frequencies,

```python
#Extract real part of the eigenvalues
Omega_squaredA = eValues.real[0]
Omega_squaredB = eValues.real[1]

# Extract the modal frequencies
OmegaA = math.sqrt(Omega_squaredA)
OmegaB = math.sqrt(Omega_squaredB)
```

And the mode shapes,

```python
# Extract the eigenvectors for each eigenvalue
eVecA = eVectors[:,0].reshape(2,1)
eVecB = eVectors[:,1].reshape(2,1)

# Assign a value of 1 to first element in each eigenvector and scale appropriately
eVecA1 = eVecA/eVecA[0]
eVecB1 = eVecB/eVecB[0]
```

Summarising our outputs,

```python
print('The first modal frequency is {one} rads/s and has a mode shape of:'
.format(one=round(OmegaB,2)))
print(eVecB1)

print('The second modal frequency is {one} rads/sec and has a mode shape of:'
.format(one=round(OmegaA,2)))
print(eVecA1)
```

- The first modal frequency is $7.0\: rads/s$ and has a mode shape of:

$$
\begin{Bmatrix}
1\\
1.45
\end{Bmatrix}
$$

- The second modal frequency is $19.46\: rads/s$ and has a mode shape of:

$$
\begin{Bmatrix}
1\\
-0.71
\end{Bmatrix}
$$

If we refer back to the last lecture we'll see that these values match the natural frequencies and mode shapes obtained by direct solution of the characteristic equation. Now that we have an easily scalable method for extracting the modal characteristics of an NDoF structure, let's apply it to the 3DoF shear building we started analysing in the previous section.

### Modal properties of our 3DoF shear building

Start by reintroducing the missing structural parameters.

```python
h3 = 4 #(m) Storey height 2-3
```

```python
m2_slab = m1_slab #(kg) Mass of level 2 concrete slab
m2_cols = h2*gamma_UC203 + h3*gamma_UC203 #(kg) Mass of half of columns above and below level 2
m2_misc = (miscAllowance/100)*(m2_slab+m2_cols) #(kg) Mass of floor beams and misc structure
m2 = m2_slab + m2_cols + m2_misc #(kg) Total lumped mass for level 2

m3_slab = m1_slab #(kg) Mass of level 3 concrete slab
m3_cols = h3*gamma_UC203 #(kg) Mass of half of columns below level 3
m3_misc = (miscAllowance/100)*(m3_slab+m3_cols) #(kg) Mass of floor beams and misc structure
m3 = m3_slab + m3_cols + m3_misc #(kg) Total lumped mass for level 3

# Mass matrix
M = np.matrix([[m1,0,0],[0, m2,0],[0,0,m3]])
```

```python
k3 = 2 * (12*E*Ix_203/h3**3) #(N/m) Combined column stiffness level 2-3

#Stiffness matrix
K = np.matrix([[k1+k2, -k2, 0],[-k2, k2+k3, -k3],[0, -k3, k3]])
```

Calculate the eigenvalues and eigenvectors,

```python
#Use library to extract eigenvalues and eigenvectors (tuple unpacking)
eValues, eVectors = la.eig(M.I*K)

#Extract real part of the eigenvalues (eVals are real numbers because matrix is symmetric)
Omega_squaredA = eValues.real[2]
Omega_squaredB = eValues.real[1]
Omega_squaredC = eValues.real[0]

# Extract the modal frequencies
OmegaA = math.sqrt(Omega_squaredA)
OmegaB = math.sqrt(Omega_squaredB)
OmegaC = math.sqrt(Omega_squaredC)

# Extract the eigenvectors for each eigenvalue
eVecA = eVectors[:,2].reshape(3,1)
eVecB = eVectors[:,1].reshape(3,1)
eVecC = eVectors[:,0].reshape(3,1)

# Assign a value of 1 to first element in each eigenvector and scale appropriately
eVecA1 = eVecA/eVecA[0]
eVecB1 = eVecB/eVecB[0]
eVecC1 = eVecC/eVecC[0]

print('The fundamental modal frequency is {one} rads/sec and has a mode shape of:'
      .format(one=round(OmegaA,2)))
print(eVecA1)

print('The second modal frequency is {one} rads/sec and has a mode shape of:'
      .format(one=round(OmegaB,2)))
print(eVecB1)

print('The third modal frequency is {one} rads/sec and has a mode shape of:'
      .format(one=round(OmegaC,2)))
print(eVecC1)
```

### Plotting the mode shapes

Again, it's worth plotting the mode shapes.

```python
#Plotting
fig, axes = plt.subplots(figsize=(15,8),nrows=1,ncols=3)
modeCoords1 = np.array([0,eVecA1[0][0], eVecA1[1][0], eVecA1[2][0]], float)
axes[0].plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-k')
axes[0].plot(modeCoords1,[0, h1, h1+h2, h1+h2+h3],'-')
axes[0].plot(modeCoords1[1],[h1],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[0].plot(modeCoords1[2],[h1+h2],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[0].plot(modeCoords1[3],[h1+h2+h3],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[0].set_xlabel('Relative amplitude')
axes[0].set_ylabel('Height (m)')
axes[0].set_title('Mode shape for $\omega_1$ =  {one} rads/sec'.format(one=round(OmegaA,1)))
axes[0].set_xlim([-2.5, 2.5])
axes[0].grid()

modeCoords2 = np.array([0,eVecB1[0][0], eVecB1[1][0], eVecB1[2][0]], float)
axes[1].plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-k')
axes[1].plot(modeCoords2,[0, h1, h1+h2, h1+h2+h3],'-')
axes[1].plot(modeCoords2[1],[h1],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[1].plot(modeCoords2[2],[h1+h2],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[1].plot(modeCoords2[3],[h1+h2+h3],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[1].set_xlabel('Relative amplitude')
axes[1].set_ylabel('Height (m)')
axes[1].set_title('Mode shape for $\omega_2$ =  {one} rads/sec'.format(one=round(OmegaB,1)))
axes[1].set_xlim([-2.5, 2.5])
axes[1].grid()

modeCoords3 = np.array([0,eVecC1[0][0], eVecC1[1][0], eVecC1[2][0]], float)
axes[2].plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-k')
axes[2].plot(modeCoords3,[0, h1, h1+h2, h1+h2+h3],'-')
axes[2].plot(modeCoords3[1],[h1],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[2].plot(modeCoords3[2],[h1+h2],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[2].plot(modeCoords3[3],[h1+h2+h3],'-o',markerfacecolor='r', markeredgecolor='k', ms=13)
axes[2].set_xlabel('Relative amplitude')
axes[2].set_ylabel('Height (m)')
axes[2].set_title('Mode shape for $\omega_3$ =  {one} rads/sec'.format(one=round(OmegaC,1)))
axes[2].set_xlim([-2.5, 2.5])
axes[2].grid()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/matrix-representation-and-the-eigenproblem/img1.jpg"
	alt="Mode shapes for mode 1, 2 and 3 (left to right). | EngineeringSkills.com"
	width="1000"
	height="685"
	caption="Fig 1. Mode shapes for mode 1, 2 and 3 (left to right)."
/>

And because we can, lets animate them too. We'll import the relevant functionality, set up a time vector and generate oscillation histories for each mass just as we did before.

```python
from matplotlib.animation import FuncAnimation
%matplotlib notebook
```

```python
# Generate a time vector
tmax = 5 #(sec) The max time
delta_t = 0.01 #(sec) The timestep
nPoints = tmax/delta_t #Number of equally spaced data points
t = np.linspace(0,tmax, int(nPoints)) # Time vector
```

```python
#Generate an oscillation for each mass in each mode

#Mode 1 oscillation
um1m1 = 1*eVecA1[0]*np.sin(OmegaA*t) #Motion of mass 1, mode 1
um2m1 = 1*eVecA1[1]*np.sin(OmegaA*t) #Motion of mass 2, mode 1
um3m1 = 1*eVecA1[2]*np.sin(OmegaA*t) #Motion of mass 3, mode 1

#Mode 2 oscillation
um1m2 = 1*eVecB1[0]*np.sin(OmegaB*t) #Motion of mass 1, mode 2
um2m2 = 1*eVecB1[1]*np.sin(OmegaB*t) #Motion of mass 2, mode 2
um3m2 = 1*eVecB1[2]*np.sin(OmegaB*t) #Motion of mass 3, mode 2

#Mode 3 oscillation
um1m3 = 1*eVecC1[0]*np.sin(OmegaC*t) #Motion of mass 1, mode 3
um2m3 = 1*eVecC1[1]*np.sin(OmegaC*t) #Motion of mass 2, mode 3
um3m3 = 1*eVecC1[2]*np.sin(OmegaC*t) #Motion of mass 3, mode 3
```

The animation code is simply repeated and extended for the additional mass.

```python
#Animation parameters
delta_t = 0.01
animLength = 10 #(sec)
frameRate = 20 #(5,10,20) frames per second (too high and animation slows down)
plotInterval = 1/frameRate #(sec) time between frame plots
dataInterval = plotInterval/delta_t #Plot moving elements every 'dataInterval-th' point

#Set up plots
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1,ncols=3, figsize=(10, 5))
ax1.set_xlabel('Relative amplitude')
ax1.set_ylabel('Height (m)')
ax1.set_title('Mode shape \n $\omega_1$ =  {one} rads/sec'.format(one=round(OmegaA,1)))
ax1.set_xlim([-2.5,2.5])

ax2.set_xlabel('Relative amplitude')
ax2.set_ylabel('Height (m)')
ax2.set_title('Mode shape \n $\omega_2$ =  {one} rads/sec'.format(one=round(OmegaB,1)))
ax2.set_xlim([-2.5,2.5])

ax3.set_xlabel('Relative amplitude')
ax3.set_ylabel('Height (m)')
ax3.set_title('Mode shape \n $\omega_3$ =  {one} rads/sec'.format(one=round(OmegaC,1)))
ax3.set_xlim([-2.5,2.5])

ax1.grid()
ax2.grid()
ax3.grid()

#Plot initial state
line1, = ax1.plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-',lw=2)
line2, = ax2.plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-',lw=2)
line3, = ax3.plot([0,0,0,0],[0, h1, h1+h2, h1+h2+h3],'-',lw=2)

dotMass1Mode1, = ax1.plot([0], h1, 'ro', markeredgecolor='k', markersize=15)
dotMass2Mode1, = ax1.plot([0], h1+h2, 'ro', markeredgecolor='k', markersize=15)
dotMass3Mode1, = ax1.plot([0], h1+h2+h3, 'ro', markeredgecolor='k', markersize=15)

dotMass1Mode2, = ax2.plot([0], h1, 'ro', markeredgecolor='k', markersize=15)
dotMass2Mode2, = ax2.plot([0], h1+h2, 'ro', markeredgecolor='k', markersize=15)
dotMass3Mode2, = ax2.plot([0], h1+h2+h3, 'ro', markeredgecolor='k', markersize=15)

dotMass1Mode3, = ax3.plot([0], h1, 'ro', markeredgecolor='k', markersize=15)
dotMass2Mode3, = ax3.plot([0], h1+h2, 'ro', markeredgecolor='k', markersize=15)
dotMass3Mode3, = ax3.plot([0], h1+h2+h3, 'ro', markeredgecolor='k', markersize=15)

#Function to animate plot objects
def animate(i):
    frm = int(i*dataInterval)

    line1.set_data([0,um1m1[frm],um2m1[frm], um3m1[frm]], [0, h1, h1+h2, h1+h2+h3])
    dotMass1Mode1.set_data(um1m1[frm],h1)
    dotMass2Mode1.set_data(um2m1[frm],h1+h2)
    dotMass3Mode1.set_data(um3m1[frm],h1+h2+h3)

    line2.set_data([0,um1m2[frm],um2m2[frm], um3m2[frm]], [0, h1, h1+h2, h1+h2+h3])
    dotMass1Mode2.set_data(um1m2[frm],h1)
    dotMass2Mode2.set_data(um2m2[frm],h1+h2)
    dotMass3Mode2.set_data(um3m2[frm],h1+h2+h3)

    line3.set_data([0,um1m3[frm],um2m3[frm], um3m3[frm]], [0, h1, h1+h2, h1+h2+h3])
    dotMass1Mode3.set_data(um1m3[frm],h1)
    dotMass2Mode3.set_data(um2m3[frm],h1+h2)
    dotMass3Mode3.set_data(um3m3[frm],h1+h2+h3)

    return

#Function to generate the animation
myAnimation = FuncAnimation(fig,
                            animate,
                            frames=int(1 + (animLength/plotInterval) ),
                            interval=plotInterval*1000,  #milliseconds
                            blit=True,
                            repeat=True)

plt.show()
myAnimation.save('3Modes.gif')
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/matrix-representation-and-the-eigenproblem/img2.gif"
	alt="Animated mode shapes for mode 1, 2 and 3 (left to right). | EngineeringSkills.com"
	width="1000"
	height="500"
	caption="Fig 2. Animated mode shapes for mode 1, 2 and 3 (left to right)."
/>


Please log in or enroll to continue

Please log in or enroll to continue