
To perform our modal analysis we can again use `SciPy` to calculate the eigenvalues and eigenvectors of the matrix $[M]^{-1}[K]$. The only additional complication is that we need to automatically sort our eigenvalues (and their corresponding eigenvectors) into ascending order. There are many ways one could achieve this and no doubt there's a more direct method than what I've done below. Feel free to investigate and implement your own method here.

Sorting the eigenvalues is pretty easy with a simple `sort()` call. To sort the eigenvectors, I started by determining the indices of the sorted eigenvalues. These are then used in a `for` loop to rearrange the columns of the modal matrix to match the order of the sorted eigenvalues. Like I said - maybe not the most efficient approach - but it works!

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

Omega = np.zeros(N) #Initialise an array of natural frequencies
phi_hat = np.matrix(np.zeros([N,N])) #Initialse an unsorted (pre-normalised) modal matrix

for i, m in enumerate(mArray):
    omega_squared =  eValues.real[i] #Extract real part of the eigenvalues
    omega = math.sqrt(omega_squared) #Extract the modal frequencies
    Omega[i] = omega #Store the natural frequency

    #Extract the eigenvectors for each eigenvalue
    eVec = eVectors[:,i].reshape(N,1)
    #Assign a value of 1 to first element in each eigenvector & scale appropriately
    eVec1 = eVec/eVec[0]
    phi_hat[:,i] = eVec1 #Add eigenvector to the modal matrix

#Sort eigenvalues and eigenvectors
sortOrder = Omega.argsort() #Sorted indices
Omega.sort() #Sort the eigenvalues (in place)

Phi_hat = np.matrix(np.zeros([N,N])) #Initialse a sorted (pre-normalised) modal matrix
for i in np.arange(N):
    Phi_hat[:,i] = phi_hat[:,sortOrder[i]]
```

### Mass normalise the modal matrix

Next is the process of mass normalisation. This becomes easy to generalise once you recognise that each column of the modal matrix is multiplied by the same scale factor. This allows us to loop through each column, calculate the scale factor for that column and apply it to the column elements in another `for` loop. Again, perhaps not the most elegant code, but it gets the job done. It's also always a good idea to confirm your normalisation has been successful by confirming $[\Phi]^T[M][\Phi] = [I]$.

```python
Phi = np.matrix(np.zeros([N,N])) #Initialse a normalised modal matrix

#For each column
for i in np.arange(N):
    #Scale factor to apply to each element in this column
    factor = np.sqrt(mArray.transpose()*np.square(Phi_hat[:,i])).item()

    #for each row element in the column
    for j in np.arange(N):
        Phi[j,i] = Phi_hat[j,i]/factor

print("[Phi].T [M] [Phi] =")
print("")
print(np.abs(np.round(Phi.transpose()*M*Phi,3)))
```

$$
[\Phi]^T[M][\Phi] =
$$

$$
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 \\
\end{bmatrix}
$$

### Visualise the vibration modes

Now we can visualise the vibration modes. Since we only have six degrees of freedom, we can generate a plot for each mode shape. However if you're analysing a structure with anything more that a dozen or so degrees of freedom, you might want to just print the first several mode shapes.

```python
#Plot size
plotHeight = 7 #Individual plot height
plotWidth = 6 #Individual plot width
```

```python
pCols = 2 if N<3 else 3 #Set number of subplot columns
pRows = int(np.ceil(N/3)) #Set number of subplot rows
fig, axes = plt.subplots(figsize=(plotWidth*plotCols,plotHeight*plotRows),
                         nrows=pRows,
                         ncols=pCols)
axes = axes.flatten()

#Plot each mode shape
for i in np.arange(N):
    xCoords = np.append(0,Phi_hat[:,i])
    xCoords = xCoords/max(abs(xCoords))
    axes[i].plot(np.zeros(N+1),floorHeights,'-k')
    axes[i].plot(xCoords,floorHeights,'-')
    axes[i].plot(xCoords[1::],floorHeights[1::],'o',
                 markerfacecolor='red',
                 markeredgecolor='k',
                 markersize=13)
    axes[i].set_ylabel('Height (m)')
    axes[i].set_title(f"$\omega$ =  {round(Omega[i],2)} rads/sec")
    axes[i].set_xlim([-1.5, 1.5])
    axes[i].grid()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/automating-modal-analysis-for-ndof/img1.jpg"
	alt="Mode shapes. | EngineeringSkills.com"
	width="1000"
	height="764"
	caption="Fig 1. Mode shapes."
/>


Please log in or enroll to continue

Please log in or enroll to continue