<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>
Frequency Response Functions (FRFs) plots
This notebook is about frequency response functions (FRFs) and the various ways they can be plotted.
Table of contents
Preamble
Dynamic system setup
Frequency response function
Nyquist plot
Bode plot
Nichols plot
Odds and ends
Preamble
We will start by setting up the computational environment for this notebook. Since it was created with Python 2.7, we will import a few things from the "future". Furthermore, we will need numpy and scipy for the numerical simulations and matplotlib for the plots
Step1:  We will also need some specific modules and a litle "IPython magic" to show the plots
Step2:  Back to top
Dynamic system setup
In this example we will simulate a two degree of freedom system (2DOF) as a LTI system. For that purpose, we will define a mass and a stiffness matrix and use proportional damping
Step3:  For the LTI system we will use a state space formulation. For that we will need the four matrices describing the system (A), the input (B), the output (C) and the feedthrough (D)
Step4:  The LTI system is simply defined as
Step5:  To check the results presented ahead we will need the angular frequencies and damping coefficients of this system. The eigenanalysis of the system matrix yields them after some computations
Step6:  Back to top
Frequency response function
A frequency response function is a complex valued function of frequency. Let us see how it looks when we plot the real and imaginary parts in separate
Step7:  Back to top
Nyquist plot
A Nyquist plot represents the real and imaginary parts of the complex FRF in a single plot
Step8:  Back to top
Bode plot
A Bode plot represents the complex FRF in magnitude-phase versus frequency
Step9:  Back to top
Nichols plot
A Nichols plot combines the Bode plot in a single plot of magnitude versus phase
<ASSISTANT_TASK:>
Python Code:
from __future__ import division, print_function
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>
Frequency Response Functions (FRFs) plots
This notebook is about frequency response functions (FRFs) and the various ways they can be plotted.
Table of contents
Preamble
Dynamic system setup
Frequency response function
Nyquist plot
Bode plot
Nichols plot
Odds and ends
Preamble
We will start by setting up the computational environment for this notebook. Since it was created with Python 2.7, we will import a few things from the "future". 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
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: We will also need some specific modules and a litle "IPython magic" to show the plots:
End of explanation
MM = np.asmatrix(np.diag([1., 2.]))
print(MM)
KK = np.asmatrix([[20., -10.],[-10., 10.]])
print(KK)
C1 = 0.1*MM+0.02*KK
print(C1)
Explanation: Back to top
Dynamic system setup
In this example we will simulate a two degree of freedom system (2DOF) as a LTI system. For that purpose, we will define a mass and a stiffness matrix and use proportional damping:
End of explanation
A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C1)]])
print(A)
Bf = KK*np.asmatrix(np.ones((2, 1)))
B = np.bmat([[np.zeros_like(Bf)],[LA.solve(MM,Bf)]])
print(B)
Cd = np.matrix((1,0))
Cv = np.asmatrix(np.zeros((1,MM.shape[1])))
Ca = np.asmatrix(np.zeros((1,MM.shape[1])))
C = np.bmat([Cd-Ca*LA.solve(MM,KK),Cv-Ca*LA.solve(MM,C1)])
print(C)
D = Ca*LA.solve(MM,Bf)
print(D)
Explanation: For the LTI system we will use a state space formulation. For that we will need the four matrices describing the system (A), the input (B), the output (C) and the feedthrough (D):
End of explanation
system = signal.lti(A, B, C, D)
Explanation: The LTI system is simply defined as:
End of explanation
w1, v1 = LA.eig(A)
ix = np.argsort(np.absolute(w1)) # order of ascending eigenvalues
w1 = w1[ix] # sorted eigenvalues
v1 = v1[:,ix] # sorted eigenvectors
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: To check the results presented ahead we will need the angular frequencies and damping coefficients of this system. The eigenanalysis of the system matrix yields them after some computations:
End of explanation
w, H = system.freqresp()
fig, ax = plt.subplots(2, 1)
fig.suptitle('Real and imaginary plots')
# Real part plot
ax[0].plot(w, H.real, label='FRF')
ax[0].axvline(wn[0], color='k', label='First mode', linestyle='--')
ax[0].axvline(wn[2], color='k', label='Second mode', linestyle='--')
ax[0].set_ylabel('Real [-]')
ax[0].grid(True)
ax[0].legend()
# Imaginary part plot
ax[1].plot(w, H.imag, label='FRF')
ax[1].axvline(wn[0], color='k', label='First mode', linestyle='--')
ax[1].axvline(wn[2], color='k', label='Second mode', linestyle='--')
ax[1].set_ylabel('Imaginary [-]')
ax[1].set_xlabel('Frequency [rad/s]')
ax[1].grid(True)
ax[1].legend()
plt.show()
Explanation: Back to top
Frequency response function
A frequency response function is a complex valued function of frequency. Let us see how it looks when we plot the real and imaginary parts in separate:
End of explanation
plt.figure()
plt.title('Nyquist plot')
plt.plot(H.real, H.imag, 'b')
plt.plot(H.real, -H.imag, 'r')
plt.xlabel('Real [-]')
plt.ylabel('Imaginary[-]')
plt.grid(True)
plt.axis('equal')
plt.show()
Explanation: Back to top
Nyquist plot
A Nyquist plot represents the real and imaginary parts of the complex FRF in a single plot:
End of explanation
w, mag, phase = system.bode()
fig, ax = plt.subplots(2, 1)
fig.suptitle('Bode plot')
# Magnitude plot
ax[0].plot(w, mag, label='FRF')
ax[0].axvline(wn[0], color='k', label='First mode', linestyle='--')
ax[0].axvline(wn[2], color='k', label='Second mode', linestyle='--')
ax[0].set_ylabel('Magnitude [dB]')
ax[0].grid(True)
ax[0].legend()
# Phase plot
ax[1].plot(w, phase*np.pi/180., label='FRF')
ax[1].axvline(wn[0], color='k', label='First mode', linestyle='--')
ax[1].axvline(wn[2], color='k', label='Second mode', linestyle='--')
ax[1].set_ylabel('Phase [rad]')
ax[1].set_xlabel('Frequency [rad/s]')
ax[1].grid(True)
ax[1].legend()
plt.show()
Explanation: Back to top
Bode plot
A Bode plot represents the complex FRF in magnitude-phase versus frequency:
End of explanation
plt.figure()
plt.title('Nichols plot')
plt.plot(phase*np.pi/180., mag)
plt.xlabel('Phase [rad/s]')
plt.ylabel('Magnitude [dB]')
plt.grid(True)
plt.show()
Explanation: Back to top
Nichols plot
A Nichols plot combines the Bode plot in a single plot of magnitude versus phase:
End of explanation
<END_TASK>