Lecture 22

E12 Linear Physical Systems Analysis

import matplotlib.pyplot as plt
from numpy import ones_like, sin, cos, linspace, logspace, sqrt, abs, log10, arctan, pi, atan2
plt.rcParams.update({'font.size': 16})

Warning

This lecture, due to technical issues, is best viewed in the ‘slides’ format available by clicking the RevealJS button on the right. You can also click here.

Frequency Response of Second Order Systems

• Second-order systems subjected to a sinusoidal input, lead to a sinusoidal output once the transients have died out.

Frequency Response of Overdamped 2nd Order Systems

Recall that overdamped means \(\zeta > 1\) where \[\zeta = \frac{b}{2\sqrt{m k}}\]

An overdamped 2nd order system

\[ \begin{aligned} m &= 50 \\ b &= 100.5 \\ k &= 1 \end{aligned} \]

• Differential equation of the system \[50 \ddot{x} + 100.5 \dot{x} + x = f(t)\]

• Transfer function: \[T(s) = \frac{X(s)}{F(s)} = \frac{1}{50s^2 + 100.5s+1}\]

• Using the roots of the denominator, we can re-write as \[T(s) = \frac{X(s)}{F(s)} = \frac{1}{(100s+1)(0.5s+1)}\]

• The (magnitude) Bode Plot needs the complex number \(\boxed{T(i\omega)}\) \[\left| T(i\omega) \right| = \left| \frac{1}{(100 i \omega +1)(0.5 i \omega + 1)}\right| = \frac{1}{\sqrt{\left( (100 \omega)^2 + 1 \right)\left( (0.5 \omega)^2 + 1 \right)}}\]

Bode Plot for overdamped 2nd order

tau1 = 100
tau2 = 0.5
omega = logspace(-4,2,500)

M_values = 1/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')
# High-frequency approximation
# ax1.loglog(omega[omega>=1/tau], 1/(omega[omega>=1/tau]*tau), ls='--',color='magenta', alpha=0.7,label="High-Frequency Approx.")
# Low-frequency approximation
# ax1.loglog(omega[omega<1/tau], ones_like(omega[omega<1/tau]), ls='--',color='blue', alpha=0.7,label="Low-Frequency Approx.")

# Corner frequencies
# ax1.axvline(1/tau1, color='red',ls=':',lw=0.5)
# ax1.axvline(1/tau2, color='red',ls=':',lw=0.5)

ax1.set_ylabel(r'$\frac{1}{\sqrt{((100\omega)^2+1)((0.5 \omega)^2+1)}}$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(False)
ax1.axhline(0, color='black')

tau1 = 100
tau2 = 0.5
omega = logspace(-4,2,500)

M_values = 1/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')
# High-frequency approximation
# ax1.loglog(omega[omega>=1/tau], 1/(omega[omega>=1/tau]*tau), ls='--',color='magenta', alpha=0.7,label="High-Frequency Approx.")
# Low-frequency approximation
# ax1.loglog(omega[omega<1/tau], ones_like(omega[omega<1/tau]), ls='--',color='blue', alpha=0.7,label="Low-Frequency Approx.")

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls=':',lw=0.5)
ax1.axvline(1/tau2, color='red',ls=':',lw=0.5)

# Label 
ax1.text(1/tau1, 1e-5, r"$\frac{1}{\tau_1} = \frac{1}{100}$", fontsize=16, color='red')
ax1.text(1/tau2, 1e-5, r"$\frac{1}{\tau_2} = \frac{1}{0.5}$", fontsize=16, color='red')


ax1.set_ylabel(r'$\frac{1}{\sqrt{((100\omega)^2+1)((0.5 \omega)^2+1)}}$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(False)
ax1.axhline(0, color='black')

• When a second-order system is overdamped, \(m s^2 + b s + k\) has two real roots.
• Each root corresponds to a corner frequency.
• We often re-write, e.g., \(50 s^2 + 100.5s + 1\) as \((\tau_1 s + 1)(\tau_2 s +1)\)

Procedure for determining \(\tau_1\) and \(\tau_2\) for overdamped systems

Suppose \(m = 2, b = 14, k = 20\). Recall that \(\zeta = \frac{b}{2\sqrt{m k}}\).

Is this system overdamped? Yes! \(\zeta = \frac{14}{2 \sqrt{2 \times 20}} \approx 1.1\)

Given the diff. eq. \(\quad 2 \ddot{x} + 14 \dot{x} + 20 x = f(t), \quad\) write the transfer function for this system in the form \(\frac{...}{(\tau_1 s + 1)(\tau_2 s + 1)}\)

\(T(s) = \frac{X}{F} = \frac{1}{2 s^2 + 14s +20}\) \(= \frac{1}{2(s+2)(s+5)}\) \(= \frac{0.05}{(0.5s+1)(0.2s+1)}\)

So the corner frequencies of this system are \(1/\tau_1 = 1/0.5 = 2\) and \(1/\tau_2 = 1/0.2 = 5\)

tau1 = 0.5
tau2 = 0.2
omega = logspace(-2,2,500)

M_values = 0.05/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')
# High-frequency approximation
# ax1.loglog(omega[omega>=1/tau], 1/(omega[omega>=1/tau]*tau), ls='--',color='magenta', alpha=0.7,label="High-Frequency Approx.")
# Low-frequency approximation
# ax1.loglog(omega[omega<1/tau], ones_like(omega[omega<1/tau]), ls='--',color='blue', alpha=0.7,label="Low-Frequency Approx.")

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1} = 2$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2} = 5$", fontsize=16, color='red')


ax1.set_ylabel(r'$\frac{0.05}{\sqrt{[(0.5\omega)^2+1][(0.2 \omega)^2+1]}}$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')

Convenient Form of Transfer Function when system is overdamped

• Starting from the differential equation \[m \ddot{x} + b \dot{x} + k x = f(t)\]
• We can write the transfer function as \[T(s) = \frac{X(s)}{F(s)} = \frac{1}{m s^2 + b s + k}\]
• and divide numerator+denominator by \(k\): \[ \frac{1/k}{(m/k) s^2 + (c/k) s + 1} = \frac{1/k}{(\tau_1 s + 1)(\tau_2 s + 1)}\]

Bode Plots for overdamped systems

Setting \(m = k = 1\) and varying only \(b\) in the range such that \(\zeta >= 1\):

# Damping ratio varied by changing b only
b = 2

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1}$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2}$", fontsize=16, color='red')

ax1.set_ylabel(r'$\left| \frac{1}{m (i \omega)^2 + b (i \omega) s + k} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

# Damping ratio varied by changing b only
b = 2.2

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1}$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2}$", fontsize=16, color='red')

ax1.set_ylabel(r'$\left| \frac{1}{m (i \omega)^2 + b (i \omega) s + k} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

# Damping ratio varied by changing b only
b = 3.0

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1}$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2}$", fontsize=16, color='red')

ax1.set_ylabel(r'$\left| \frac{1}{m (i \omega)^2 + b (i \omega) s + k} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

# Damping ratio varied by changing b only
b = 4.0

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1}$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2}$", fontsize=16, color='red')

ax1.set_ylabel(r'$\left| \frac{1}{m (i \omega)^2 + b (i \omega) s + k} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

# Damping ratio varied by changing b only
b = 10.0

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))

# --- TOP SUBPLOT ---
# Bode plots
ax1.loglog(omega, M_values,lw=2,color='black')

# Corner frequencies
ax1.axvline(1/tau1, color='red',ls='-',lw=1)
ax1.axvline(1/tau2, color='red',ls='-',lw=1)

# Label 
ax1.text(0.5/tau1, 1e-3, r"$\frac{1}{\tau_1}$", fontsize=16, color='red')
ax1.text(0.5/tau2, 1e-4, r"$\frac{1}{\tau_2}$", fontsize=16, color='red')

ax1.set_ylabel(r'$\left| \frac{1}{m (i \omega)^2 + b (i \omega) s + k} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

Dimensionless Transfer Function for 2nd order systems

\(\displaystyle T(s) = \frac{X(s)}{F(s)} = \frac{1}{m s^2 + b s + k}\)

\(\displaystyle T(s) = \frac{k X(s)}{F(s)} = \frac{k }{m s^2 + b s + k}\) \(= \frac{k/k}{(m/k) s^2 + (b/k) s + (k/k)}\)

\(\displaystyle T(s) = \frac{1}{(m/k) s^2 + (b/k) s + 1} = \frac{1}{\left(\frac{s}{\omega_n}\right)^2 + 2 \zeta (s/\omega_n) + 1}\)

\(\displaystyle T(s) = \frac{k X(s)}{F(s)} = \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\)

• This transfer function is unitless
• It is a function of two quantities (\(\zeta, \omega_n\)) instead of three (\(m,b,k\))

Bode Plot of general 2nd-order systems (overdamped and underdamped)

• Starting from the transfer function \(T(s) = \frac{\omega_n^2}{s^2+2 \zeta \omega_n s + \omega_n^2}\)

\(T(i\omega) = \frac{\omega_n^2}{(i \omega)^2 + 2 \zeta \omega_n (i \omega) + \omega_n^2}\) \(\phantom{T(i\omega) }= \frac{\omega_n^2}{- \omega^2 + \omega_n^2 + 2 \zeta \omega_n \omega i}\)

• We introduce a new quantity \(\boxed{r = \frac{\omega}{\omega_n}}\): the ratio of the input frequency \(\omega\) to the undamped natural frequency \(\omega_n\).

\(T(i\omega) = \frac{\omega_n^2}{- \left( \frac{\omega}{\omega_n} \right)^2 \cdot \omega_n^2 + \omega_n^2 + 2 \zeta \omega_n \left( \frac{\omega}{\omega_n} \right) \cdot \omega_n i }\)

\(\phantom{T(i\omega)} = \frac{\omega_n^2}{\omega_n^2 \left(1 - r^2 \right) + 2 \zeta r \omega_n^2 i}\) \(= \frac{1}{\left(1 - r^2 \right) + 2 \zeta r i}\)

The magnitude of this complex number is \(\frac{1}{\sqrt{(1-r)^2 + 4 \zeta^2 r^2}}\)

Bode plot as a function of frequency ratio

The (magnitude) Bode plot for a 2nd order system \(T(s) = \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\) can therefore be represented as \[|T(i \omega)| = \frac{1}{\sqrt{(1-r^2)^2 + 4 \zeta^2 r^2}}\]

# Damping ratio varied by changing b only
b = 10.0

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(omega, M_values,lw=2,color='black')

ax1.set_ylabel(r'Gain $\left| \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2} \right|$')
ax1.set_xlabel('Frequency of input $\omega$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\omega_n = 1, \quad \zeta = {b/2}$")
plt.show()

# Damping ratio varied by changing b only
b = 10.0

# m and k set to 1 for simplicity
k = 1
m = 1

# Formula for roots of polynomial when k = m = 1
s1 = 0.5*(-b-sqrt(b**2-4)) 
s2 = 0.5*(-b+sqrt(b**2-4))
tau1 = -s1
tau2 = -s2

omega = logspace(-4,2,500)

M_values = (1/k)/sqrt((1+(omega*tau1)**2)*(1+(omega*tau2)**2))
fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(omega, M_values,lw=2,color='black')

ax1.set_ylabel(r'Gain $\frac{1}{\sqrt{(1-r^2)^2 + 4 \zeta^2 r^2}}$')
ax1.set_xlabel('Ratio of input frequency to $\omega_n$,  $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.set_title(f"$\quad \zeta = {b/2}$")
plt.show()

• In this form, the Bode plot has \(r\) on the horizontal axis, gain on the vertical axis, and a single parameter \(\zeta\).

Plotting \(T(i \omega)\) for 2nd order systems

Examples of 2nd order Bode Plots

# Damping ratio varied by changing b only
z = 2

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 1.0

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 0.8

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 0.6

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 0.3

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 0.1

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z = 0.01

r = logspace(-2,2,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

Complex Zeros of the Transfer function

• If a system has complex zeros, then its Bode plot looks similar to the case when it has complex poles, except that it is flipped upside down.

tau1 = 5
omega = logspace(-2,2,500)

M_values1 = sqrt(9 - 5*omega**2 + omega**4)/3

fig, ax1 = plt.subplots(1, 1, figsize=(6, 3))

# --- TOP SUBPLOT ---
ax1.loglog(omega, M_values1, label=f"$T(s) = (s^2 + s + 3)/3$")
ax1.set_title(f'A system with a pair of complex zeros')
ax1.set_ylabel(f"Gain $|T(i\omega)|$")
ax1.set_xlabel('Frequency of input $\\omega$')
ax1.grid(True, which="both", ls="-", alpha=0.5)
ax1.axhline(0, color='black')
ax1.scatter(1/tau1,sqrt(2))
plt.tight_layout()
plt.legend(loc="upper left",fontsize=16)

What happens at v. low \(\zeta\)

Let’s zoom in to the Bode plot at small damping ratios.

# Damping ratio varied by changing b only
z = 0.01

r = logspace(-0.5,0.5,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-1,1e5])
plt.show()

# Damping ratio varied by changing b only
z = 0.001

r = logspace(-0.5,0.5,500)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-1,1e5])
plt.show()

# Damping ratio varied by changing b only
z = 0.0001

r = logspace(-0.5,0.5,2000)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-1,1e5])
plt.show()

# Damping ratio varied by changing b only
z = 0

r = logspace(-0.5,0.5,2000)

M_values = 1/sqrt((1-r**2)**2 + 4*(z*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values,lw=2,color='black',label=f"$\zeta = {z}$")
ax1.axvline(x=1,ymin=0.65,ymax=1, color='black',ls='--',lw=2)
ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-1,1e5])
plt.show()

Resonance

• We notice that there is some critical value of \(\zeta\) below which the Bode plot has a ‘hump’.
• This phenomenon is known as resonance: when the system \(T(s) = \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\) has Gain \(>1\).

# Damping ratio varied by changing b only
z1 = 1.2

r = logspace(-2,2,500)

M_values1 = 1/sqrt((1-r**2)**2 + 4*(z1*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values1,lw=2,color='black',label=f"$\zeta = {z1}$")


ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z1 = 1.2
z2 = 0.8

r = logspace(-2,2,500)

M_values1 = 1/sqrt((1-r**2)**2 + 4*(z1*r)**2)
M_values2 = 1/sqrt((1-r**2)**2 + 4*(z2*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values1,lw=2,color='black',label=f"$\zeta = {z1}$")
ax1.loglog(r, M_values2,lw=2,color='gray',label=f"$\zeta = {z2}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

# Damping ratio varied by changing b only
z1 = 1.2
z2 = 0.8
z3 = 0.6

r = logspace(-2,2,500)

M_values1 = 1/sqrt((1-r**2)**2 + 4*(z1*r)**2)
M_values2 = 1/sqrt((1-r**2)**2 + 4*(z2*r)**2)
M_values3 = 1/sqrt((1-r**2)**2 + 4*(z3*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values1,lw=2,color='black',label=f"$\zeta = {z1}$")
ax1.loglog(r, M_values2,lw=2,color='gray',label=f"$\zeta = {z2}$")
ax1.loglog(r, M_values3,lw=2,color='red',label=f"$\zeta = {z3}$")

ax1.set_ylabel(r'Gain')
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

When does resonance happen?

# Damping ratio varied by changing b only
z1 = 0.4
z2 = 0.2
z3 = 0.05

r = logspace(-2,2,500)

M_values1 = 1/sqrt((1-r**2)**2 + 4*(z1*r)**2)
M_values2 = 1/sqrt((1-r**2)**2 + 4*(z2*r)**2)
M_values3 = 1/sqrt((1-r**2)**2 + 4*(z3*r)**2)

fig, ax1 = plt.subplots(1,1, figsize=(8, 4))
ax1.loglog(r, M_values1,lw=2,color='xkcd:beige',label=f"$\zeta = {z1}$")
ax1.loglog(r, M_values2,lw=2,color='xkcd:salmon',label=f"$\zeta = {z2}$")
ax1.loglog(r, M_values3,lw=2,color='xkcd:sky blue',label=f"$\zeta = {z3}$")

ax1.set_ylabel(r'Gain')
ax1.axvline(1, color='black',ls='--',lw=1.0)
ax1.set_xlabel('Frequency ratio $r = \omega/\omega_n$')
ax1.grid(True,which='minor',ls=':',lw=0.5)
ax1.axhline(0, color='black')
ax1.legend()
ax1.set_title("Magnitude Bode Plot")
plt.ylim([1e-2,1e1])
plt.show()

• A maximum in \(\frac{1}{\sqrt{(1-r^2)^2 + 4 \zeta^2 r^2}}\) only occurs when \(r = \sqrt{1 - 2\zeta^2}\)
• If \(\zeta > \frac{\sqrt{2}}{2} \approx 0.707\), resonance does not occur.
• For \(0 < \zeta < \frac{\sqrt{2}}{2}\), the resonant frequency is \[\omega_r = \omega_n \sqrt{1-2\zeta^2}\]
  • If the system is excited at this frequency, it will experience resonance.
• The resonant gain \(T(i \omega_r)\) can be shown to be \[|T(i \omega_r)| = \frac{1}{2 \zeta \sqrt{1- \zeta^2}}\]

Comprehensive Bode Plot for second-order systems

z1 = 0.01
z2 = 0.1
z3 = 0.5
z4 = 0.7
z5 = 1.0

r = logspace(-2,2,500)

M_values1 = 1/sqrt((1-r**2)**2 + 4*(z1*r)**2)
M_values2 = 1/sqrt((1-r**2)**2 + 4*(z2*r)**2)
M_values3 = 1/sqrt((1-r**2)**2 + 4*(z3*r)**2)
M_values4 = 1/sqrt((1-r**2)**2 + 4*(z4*r)**2)
M_values5 = 1/sqrt((1-r**2)**2 + 4*(z5*r)**2)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 9))

phi_values1 = -atan2((2*z1*r),(1-r**2))
phi_values2 = -atan2((2*z2*r),(1-r**2))
phi_values3 = -atan2((2*z3*r),(1-r**2))
phi_values4 = -atan2((2*z4*r),(1-r**2))
phi_values5 = -atan2((2*z5*r),(1-r**2))

# --- TOP SUBPLOT ---

ax1.loglog(r, M_values1, label=f"$\\zeta = {z1}$")
ax1.loglog(r, M_values2, label=f"$\\zeta = {z2}$")
ax1.loglog(r, M_values3, label=f"$\\zeta = {z3}$")
ax1.loglog(r, M_values4, label=f"$\\zeta = {z4}$")
ax1.loglog(r, M_values5, label=f"$\\zeta = {z5}$")
ax1.set_title('Second Order Systems')
ax1.set_ylabel('Gain ')
ax1.grid(True, which="both", ls="-", alpha=0.5)
ax1.legend(loc="lower left")
ax1.axhline(0, color='black')
#ax1.axvline(1/tau, color='red',ls='--')
#ax1.axvline(1/tau2, color='red',ls='--')

# --- BOTTOM SUBPLOT (Identical) ---
ax2.semilogx(r, phi_values1, label=f"$\\zeta = {z1}$")
ax2.semilogx(r, phi_values2, label=f"$\\zeta = {z2}$")
ax2.semilogx(r, phi_values3, label=f"$\\zeta = {z3}$")
ax2.semilogx(r, phi_values4, label=f"$\\zeta = {z4}$")
ax2.semilogx(r, phi_values5, label=f"$\\zeta = {z5}$")
ax2.set_xlabel('Frequency Ratio $r = \omega / \omega_n$')
ax2.set_ylabel('Phase Shift ')
ticks = [0, -pi/6, -pi/3, -pi/2, -2*pi/3, -5*pi/6, -pi]

labels = [
    '0', 
    r'$-30^{\circ}$', 
    r'$-60^{\circ}$', 
    r'$-90^{\circ}$', 
    r'$-120^{\circ}$', 
    r'$-150^{\circ}$', 
    r'$-180^{\circ}$'
]
ax2.set_yticks(ticks)
ax2.set_yticklabels(labels)
ax2.grid(True, which="both", ls="-", alpha=0.5)
ax2.legend(loc="lower left")
ax2.axhline(0, color='black')

plt.tight_layout()

Resonant frequency vs. natural frequency

Three frequencies inherent to a second-order system.

The undamped natural frequency \[\omega_n = \sqrt{k/m}\] at which the system would hypothetically oscillate if there were no damping.

The damped natural frequency \[\omega_d = \omega_n \sqrt{1-\zeta^2}\] at which an underdamped system actually oscillates.

The resonant frequency \[\omega_r = \omega_n \sqrt{1-2\zeta^2}\] at which, if we provide an input to the system, the system’s oscillations will grow with time.

In-class activity

Predict the frequency with which the spring-mass system from Tuesday’s class should be oscillated to achieve resonance. Give your answer in Hertz.

Note that each mass was 200 grams, and the spring is listed as having a spring constant of \(0.2~\mathrm{lbf}/\mathrm{in}\).