For an LTI system with transfer function \( G(s) \): a sinusoidal input \( u(t)=U\cos(\omega t) \): produces a steady-state output \( y(t)=U|G(j\omega)|\cos(\omega t+\angle G(j\omega)) \):. The map \( \omega \mapsto G(j\omega) \): is the frequency response.

Resonance occurs when \( |G(j\omega)| \): has a local maximum for some \( \omega > 0 \):. Bandwidth quantifies how rapidly the response rolls off around a key frequency (often around the resonance), and the quality factor \( Q \): quantifies the sharpness of resonance (high \( Q \): means a narrow and tall resonant peak).

Throughout, we use angular frequency \( \omega \): in rad/s. The conversion to Hertz is \( f = \omega/(2\pi) \):.

flowchart TD
  A["Physical system (mechanical or electrical)"] --> B["Linear model (ODE)"]
  B --> C["Transfer function G(s)"]
  C --> D["Evaluate on j*w : G(j*w)"]
  D --> E["Magnitude |G(j*w)| and phase angle(G(j*w))"]
  E --> F["Resonance peak (if any)"]
  F --> G["Half-power points (w1,w2)"]
  G --> H["Bandwidth BW = w2 - w1 and quality factor Q = w0/BW"]
  H --> I["Interpretation: damping, energy loss, design trade-offs"]
        

A large class of mechanical and electrical subsystems reduce (after normalization) to a second-order transfer function. We will use the normalized low-pass form (unit DC gain):

\[ G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad \omega_n > 0, \quad \zeta \ge 0. \]

Let \( r = \omega/\omega_n \):. By direct substitution \( s=j\omega \): and using \( |a+jb|^2=a^2+b^2 \):, the magnitude-squared is

\[ |G(j\omega)|^2 = \frac{1}{\left(1-r^2\right)^2 + \left(2\zeta r\right)^2}. \]

Proposition 1 (Existence and location of resonant peak). The function \( |G(j\omega)| \): has an interior maximum at \( \omega_r = \omega_n\sqrt{1-2\zeta^2} \): if and only if \( \zeta < 1/\sqrt{2} \):.

Proof. Maximizing \( |G(j\omega)|^2 \): is equivalent to minimizing the denominator \( D(r) = (1-r^2)^2 + (2\zeta r)^2 \):. Expand:

\[ D(r) = (1-2r^2+r^4) + 4\zeta^2 r^2 = r^4 + (4\zeta^2-2)r^2 + 1. \]

Differentiate w.r.t. \( r \):

\[ \frac{dD}{dr} = 4r^3 + 2(4\zeta^2-2)r = 4r\left(r^2 - (1-2\zeta^2)\right). \]

Critical points satisfy \( r=0 \): or \( r^2 = 1-2\zeta^2 \):. A positive interior critical point exists iff \( 1-2\zeta^2 > 0 \):, i.e. \( \zeta < 1/\sqrt{2} \):, giving \( r_r = \sqrt{1-2\zeta^2} \): and \( \omega_r = \omega_n r_r \):. Since \( D(r) \to \infty \): as \( r \to \infty \):, this interior point is the unique global minimum of \( D(r) \): on \( (0,\infty) \): and hence the unique resonant peak. ∎

Evaluating the denominator at \( r_r^2 = 1-2\zeta^2 \): yields the peak magnitude:

\[ |G(j\omega_r)| = \frac{1}{2\zeta\sqrt{1-\zeta^2}}, \quad \text{valid when } \zeta < 1/\sqrt{2}. \]

When \( \zeta \ge 1/\sqrt{2} \):, no interior peak exists and the magnitude decreases monotonically from its DC value \( |G(0)|=1 \):.

For resonant responses, a standard quantitative definition uses half-power frequencies \( \omega_1 < \omega_r < \omega_2 \): defined by the condition \( |G(j\omega)|^2 = \tfrac{1}{2}|G(j\omega_r)|^2 \):, i.e. a −3 dB drop relative to the peak. Define the bandwidth as \( \mathrm{BW} = \omega_2 - \omega_1 \): and the quality factor as \( Q = \omega_0/\mathrm{BW} \): where \( \omega_0 \): is a characteristic frequency (often \( \omega_0=\omega_n \): or \( \omega_0=\omega_r \): for light damping).

Proposition 2 (Exact half-power frequencies for normalized second order). Assume \( \zeta < 1/\sqrt{2} \):. Let \( r = \omega/\omega_n \):. The half-power points relative to the peak satisfy

\[ r_{1,2}^2 = 1 - 2\zeta^2 \mp 2\zeta\sqrt{1-\zeta^2}, \quad \omega_{1,2} = \omega_n\sqrt{r_{1,2}^2}. \]

Proof sketch. The peak occurs at \( r_r^2 = 1-2\zeta^2 \):, and the minimum denominator value is

\[ D_{\min} = D(r_r) = 4\zeta^2(1-\zeta^2). \]

The half-power condition \( |G(j\omega)|^2 = \tfrac{1}{2}|G(j\omega_r)|^2 \): is equivalent to \( D(r) = 2D_{\min} \):. Writing \( y=r^2 \):, the equation becomes a quadratic whose discriminant simplifies to \( 16\zeta^2(1-\zeta^2) \):, producing the stated roots. ∎

Corollary 1 (Lightly damped approximation). If \( 0 < \zeta \ll 1 \): then \( \omega_1 \approx \omega_n(1-\zeta) \):, \( \omega_2 \approx \omega_n(1+\zeta) \):, so

\[ \mathrm{BW} = \omega_2-\omega_1 \approx 2\zeta\omega_n, \quad Q \approx \frac{\omega_n}{\mathrm{BW}} \approx \frac{1}{2\zeta}. \]

This approximation is the main reason that \( Q \): is often treated as the inverse of damping.

Remark (Control bandwidth vs resonant bandwidth). In control engineering, “bandwidth” can also mean the frequency at which a low-pass closed-loop magnitude first drops to \( 1/\sqrt{2} \): of its low-frequency value (−3 dB point). In this lesson, “bandwidth” refers to the resonant half-power bandwidth around the resonant peak, which is the quantity that connects most directly to \( Q \):.

A physically transparent definition of quality factor is

\[ Q = 2\pi\,\frac{E_{\text{stored}}}{\Delta E_{\text{dissipated per cycle}}}, \quad \text{evaluated near resonance.} \]

Consider the homogeneous underdamped oscillator \( x'' + 2\zeta\omega_n x' + \omega_n^2 x = 0 \): with \( 0 < \zeta < 1 \):. Its solution is \( x(t)=A e^{-\zeta\omega_n t}\cos(\omega_d t-\phi) \): where \( \omega_d = \omega_n\sqrt{1-\zeta^2} \):.

Define total energy (per unit mass, up to scaling) as \( E(t) = \tfrac{1}{2}x'(t)^2 + \tfrac{1}{2}\omega_n^2 x(t)^2 \):. A standard averaging argument over one period shows the envelope decays as

\[ E(t) \approx E(0) e^{-2\zeta\omega_n t}. \]

Over one cycle \( T = 2\pi/\omega_d \):, the fractional energy loss is

\[ \frac{\Delta E}{E} = 1 - e^{-2\zeta\omega_n T} = 1 - e^{-4\pi\zeta\,\omega_n/\omega_d}. \]

For light damping \( \zeta \ll 1 \):, use \( e^{-\epsilon} \approx 1-\epsilon \): and \( \omega_n/\omega_d \approx 1 \): to obtain

\[ \Delta E \approx E\,4\pi\zeta \quad \Longrightarrow \quad Q = 2\pi\,\frac{E}{\Delta E} \approx \frac{1}{2\zeta}. \]

Therefore, the energy definition of \( Q \): is consistent with both the half-power bandwidth definition and the pole-based interpretation of damping.

For forced vibration \( m\ddot x + c\dot x + kx = f(t) \):, the force-to-displacement transfer function is

\[ \frac{X(s)}{F(s)} = \frac{1}{ms^2 + cs + k} = \frac{1}{k}\,\frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad \omega_n = \sqrt{\frac{k}{m}},\; \zeta = \frac{c}{2\sqrt{km}}. \]

Hence the resonance condition and resonant frequency are inherited from the normalized model, but the overall magnitude is scaled by \( 1/k \):. In particular, when \( \zeta < 1/\sqrt{2} \):,

\[ \left|\frac{X(j\omega_r)}{F(j\omega_r)}\right| = \frac{1}{k}\,\frac{1}{2\zeta\sqrt{1-\zeta^2}}. \]

Interpretation. Small viscous damping \( c \): yields small \( \zeta \):, hence large \( Q \): and a sharp resonance. In mechanical design this implies large vibration amplitudes under narrow-band excitation; in control design it implies strong sensitivity to forcing at particular frequencies and motivates damping augmentation or notch filtering (covered later when shaping frequency response).

RLC circuits exhibit the same second-order structure. Let \( \omega_0 = 1/\sqrt{LC} \):.

Series RLC (band-pass current resonance). For a series connection, the impedance is \( Z(j\omega)=R + j(\omega L - 1/(\omega C)) \):. The current magnitude is

\[ |I(j\omega)| = \frac{|V(j\omega)|}{\sqrt{R^2 + \left(\omega L - \frac{1}{\omega C}\right)^2}}. \]

The peak occurs at \( \omega_0 \):. The half-power condition relative to the peak \( |I| = |V|/(\sqrt{2}R) \): yields \( \left|\omega L - \tfrac{1}{\omega C}\right| = R \):. Multiplying by \( \omega \): and solving the two quadratics gives

\[ \omega_{1} = \frac{-R + \sqrt{R^2 + 4\frac{L}{C}}}{2L}, \quad \omega_{2} = \frac{R + \sqrt{R^2 + 4\frac{L}{C}}}{2L}, \quad \mathrm{BW} = \omega_2 - \omega_1 = \frac{R}{L}. \]

Therefore the series-RLC quality factor is exactly

\[ Q = \frac{\omega_0}{\mathrm{BW}} = \frac{\omega_0 L}{R} = \frac{1}{\omega_0 C R}. \]

Parallel RLC (impedance resonance). For a parallel connection, the admittance is \( Y(j\omega)=\tfrac{1}{R} + j(\omega C - 1/(\omega L)) \):. The impedance magnitude peaks at \( \omega_0 \): with \( |Z(j\omega_0)|=R \):. The half-power condition yields \( |\omega C - 1/(\omega L)| = 1/R \):, leading to

\[ \mathrm{BW} = \omega_2 - \omega_1 = \frac{1}{RC}, \quad Q = \frac{\omega_0}{\mathrm{BW}} = \omega_0 RC = \frac{R}{\omega_0 L}. \]

These exact relationships illustrate why \( Q \): is both a frequency-domain metric (sharpness of resonance) and an energy metric (loss in the resistor per cycle).

In practice, you often have sampled frequency-response magnitude data (from simulation or experiments). For a system that is well-approximated by a dominant second-order resonance, you can estimate \( \omega_r \): from the peak location, then estimate bandwidth from half-power points, then infer \( Q \): and (for lightly damped cases) an estimate of \( \zeta \approx 1/(2Q) \):.

flowchart TD
  A["Input: sampled magnitude data {w_i, |G(j*w_i)|}"] --> B["Find peak index i* (max magnitude)"]
  B --> C["wr_hat = w_i* and Mr_hat = |G(j*wr_hat)|"]
  C --> D["Target = Mr_hat / sqrt(2)"]
  D --> E["Search left for first crossing => w1_hat"]
  D --> F["Search right for first crossing => w2_hat"]
  E --> G["BW_hat = w2_hat - w1_hat"]
  F --> G
  G --> H["Q_hat = wr_hat / BW_hat"]
  H --> I["If lightly damped: zeta_hat approx 1/(2*Q_hat)"]
        

When the resonance is not isolated (multiple modes, zeros near the peak, strong sensor/actuator dynamics), the correct approach is frequency-domain model fitting (e.g., least-squares on complex data) or state-space identification, introduced in later chapters.

The implementations below (i) compute \( \omega_n, \zeta, Q, \omega_r \): analytically for canonical models, (ii) compute half-power bandwidth (exact and approximate), and (iii) estimate resonant parameters from sampled frequency-response data.

Chapter12_Lesson3.py

Python libraries commonly used in system dynamics include numpy, scipy.signal, and the control-oriented package control (python-control) for classical frequency-response workflows.

""" 
Chapter12_Lesson3.py
System Dynamics — Chapter 12, Lesson 3
Resonance, Bandwidth, and Quality Factor in Mechanical and Electrical Systems

Dependencies (recommended):
  - numpy
  - scipy
Optional:
  - control  (python-control)
"""

from __future__ import annotations
import math
from dataclasses import dataclass
from typing import Optional, Tuple, Sequence

import numpy as np

try:
    from scipy import signal
except Exception as e:
    signal = None


@dataclass(frozen=True)
class SecondOrderParams:
    wn: float     # natural frequency (rad/s)
    zeta: float   # damping ratio
    Q: float      # quality factor (dimensionless)


def second_order_from_mck(m: float, c: float, k: float) -> SecondOrderParams:
    """
    Mass-spring-damper: m x¨ + c x˙ + k x = f
    Transfer from force to displacement: X/F = 1 / (m s^2 + c s + k)

    wn = sqrt(k/m)
    zeta = c / (2*sqrt(k*m))
    Q = 1/(2*zeta)   (valid for standard 2nd-order form)
    """
    if m <= 0 or k <= 0:
        raise ValueError("m and k must be positive.")
    wn = math.sqrt(k / m)
    zeta = c / (2.0 * math.sqrt(k * m))
    if zeta <= 0:
        Q = float("inf")
    else:
        Q = 1.0 / (2.0 * zeta)
    return SecondOrderParams(wn=wn, zeta=zeta, Q=Q)


def resonance_frequency(wn: float, zeta: float) -> Optional[float]:
    """
    Resonant frequency (magnitude peak) for the normalized low-pass 2nd-order:
        G(s) = wn^2 / (s^2 + 2 zeta wn s + wn^2)
    exists only if zeta < 1/sqrt(2)
        wr = wn * sqrt(1 - 2 zeta^2)
    """
    if zeta < 0:
        raise ValueError("zeta must be nonnegative.")
    if zeta >= 1.0 / math.sqrt(2.0):
        return None
    return wn * math.sqrt(1.0 - 2.0 * zeta * zeta)


def peak_magnitude_normalized_lowpass(zeta: float) -> Optional[float]:
    """
    Peak |G(jw)| for normalized low-pass second-order above.
    If zeta >= 1/sqrt(2), no interior peak; max occurs at w=0 and equals 1.
    If zeta < 1/sqrt(2):
        Mr = 1 / (2 zeta sqrt(1 - zeta^2))
    """
    if zeta >= 1.0 / math.sqrt(2.0):
        return None
    return 1.0 / (2.0 * zeta * math.sqrt(1.0 - zeta * zeta))


def half_power_frequencies_relative_to_peak(wn: float, zeta: float) -> Optional[Tuple[float, float]]:
    """
    Half-power frequencies w1 < w2 around the resonant peak of normalized low-pass.
    Defined by: |G(jw)|^2 = (1/2) |G(jwr)|^2

    Exact (in terms of r = w/wn):
        r1^2 = 1 - 2 zeta^2 - 2 zeta sqrt(1 - zeta^2)
        r2^2 = 1 - 2 zeta^2 + 2 zeta sqrt(1 - zeta^2)

    Requires zeta < 1/sqrt(2) (peak exists).
    """
    wr = resonance_frequency(wn, zeta)
    if wr is None:
        return None
    inside = 1.0 - zeta * zeta
    if inside <= 0:
        return None
    r1_sq = 1.0 - 2.0 * zeta * zeta - 2.0 * zeta * math.sqrt(inside)
    r2_sq = 1.0 - 2.0 * zeta * zeta + 2.0 * zeta * math.sqrt(inside)
    if r1_sq <= 0 or r2_sq <= 0:
        return None
    return (wn * math.sqrt(r1_sq), wn * math.sqrt(r2_sq))


def bandwidth_from_half_power(w1: float, w2: float) -> float:
    return float(w2 - w1)


def approximate_bandwidth_small_damping(wn: float, zeta: float) -> float:
    """
    For zeta << 1:
        Delta w ≈ 2 zeta wn
    """
    return 2.0 * zeta * wn


def estimate_peak_and_half_power_from_samples(w: Sequence[float], mag: Sequence[float]) -> dict:
    """
    Estimate resonant peak and half-power points from sampled magnitude data.

    Inputs:
      w   : frequencies in rad/s (monotone increasing)
      mag : |G(jw)| magnitudes at corresponding points

    Outputs:
      dict with keys: wr_hat, Mr_hat, w1_hat, w2_hat, bw_hat, Q_hat (if possible)
    """
    w = np.asarray(w, dtype=float)
    mag = np.asarray(mag, dtype=float)
    if w.ndim != 1 or mag.ndim != 1 or w.size != mag.size:
        raise ValueError("w and mag must be 1D arrays of same length.")
    if np.any(np.diff(w) <= 0):
        raise ValueError("w must be strictly increasing.")

    idx = int(np.argmax(mag))
    wr_hat = float(w[idx])
    Mr_hat = float(mag[idx])

    target = Mr_hat / math.sqrt(2.0)  # half-power magnitude relative to peak

    # Find left crossing
    w1_hat = None
    for i in range(idx, 0, -1):
        if (mag[i] - target) * (mag[i-1] - target) <= 0:
            # linear interpolation
            t = (target - mag[i-1]) / (mag[i] - mag[i-1] + 1e-300)
            w1_hat = float(w[i-1] + t * (w[i] - w[i-1]))
            break

    # Find right crossing
    w2_hat = None
    for i in range(idx, w.size - 1):
        if (mag[i] - target) * (mag[i+1] - target) <= 0:
            t = (target - mag[i]) / (mag[i+1] - mag[i] + 1e-300)
            w2_hat = float(w[i] + t * (w[i+1] - w[i]))
            break

    out = {"wr_hat": wr_hat, "Mr_hat": Mr_hat, "w1_hat": w1_hat, "w2_hat": w2_hat}
    if w1_hat is not None and w2_hat is not None and w2_hat > w1_hat:
        bw = w2_hat - w1_hat
        out["bw_hat"] = float(bw)
        out["Q_hat"] = float(wr_hat / bw) if bw > 0 else None
    else:
        out["bw_hat"] = None
        out["Q_hat"] = None
    return out


def demo_mass_spring_damper():
    # Example physical parameters (units: kg, N*s/m, N/m)
    m, c, k = 1.0, 0.4, 100.0
    p = second_order_from_mck(m, c, k)

    wr = resonance_frequency(p.wn, p.zeta)
    hw = half_power_frequencies_relative_to_peak(p.wn, p.zeta)
    print("Mass–spring–damper parameters:")
    print(f"  wn   = {p.wn:.6g} rad/s")
    print(f"  zeta = {p.zeta:.6g}")
    print(f"  Q    = {p.Q:.6g}")
    print(f"  wr   = {wr:.6g} rad/s" if wr else "  wr   = (no resonant peak)")
    if hw:
        w1, w2 = hw
        bw = bandwidth_from_half_power(w1, w2)
        print(f"  w1   = {w1:.6g} rad/s")
        print(f"  w2   = {w2:.6g} rad/s")
        print(f"  BW   = {bw:.6g} rad/s")
        print(f"  Q_hp = {wr/bw:.6g} (wr/BW)")
        print(f"  BW approx (2 zeta wn) = {approximate_bandwidth_small_damping(p.wn, p.zeta):.6g} rad/s")

    # Numerical frequency response check
    if signal is None:
        print("\nscipy is not available; skipping numerical frequency response demo.")
        return

    # Force -> displacement transfer: G(s) = 1/(m s^2 + c s + k)
    sys = signal.TransferFunction([1.0], [m, c, k])
    w = np.logspace(-1, 3, 4000)  # rad/s
    w, H = signal.freqresp(sys, w=w)
    mag = np.abs(H)

    est = estimate_peak_and_half_power_from_samples(w, mag)
    print("\nEstimated from sampled frequency response:")
    for k, v in est.items():
        print(f"  {k}: {v}")


def demo_series_rlc():
    """
    Series RLC band-pass across R:
      Z = R + j(wL - 1/(wC))
      I = V / Z
      Vr = I*R -> |Vr/V| = R/|Z|

    Resonance at w0 = 1/sqrt(LC)
    Quality factor (series): Q = w0 L / R = 1/(w0 C R)
    Bandwidth (rad/s): Delta w = R/L
    """
    R, L, C = 10.0, 50e-3, 10e-6
    w0 = 1.0 / math.sqrt(L * C)
    Q = w0 * L / R
    bw = R / L
    print("\nSeries RLC example:")
    print(f"  w0 = {w0:.6g} rad/s")
    print(f"  Q  = {Q:.6g}")
    print(f"  BW = {bw:.6g} rad/s")


if __name__ == "__main__":
    demo_mass_spring_damper()
    demo_series_rlc()
      

Chapter12_Lesson3.cpp

In C++, a typical workflow is to use std::complex for frequency-response evaluation and pair it with Eigen (linear algebra) or Boost (numerical utilities) as systems scale.

// Chapter12_Lesson3.cpp
// System Dynamics — Chapter 12, Lesson 3
// Resonance, Bandwidth, and Quality Factor in Mechanical and Electrical Systems
//
// Build (example):
//   g++ -O2 -std=c++17 Chapter12_Lesson3.cpp -o Chapter12_Lesson3
//
// Notes:
// - Uses only the C++ standard library (std::complex) for transfer-function evaluation.
// - For larger projects, pair this with Eigen (linear algebra) and/or Boost (numerics).

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
#include <limits>

struct SecondOrderParams {
    double wn;   // rad/s
    double zeta; // damping ratio
    double Q;    // quality factor
};

SecondOrderParams second_order_from_mck(double m, double c, double k) {
    if (m <= 0.0 || k <= 0.0) throw std::runtime_error("m and k must be positive.");
    double wn = std::sqrt(k / m);
    double zeta = c / (2.0 * std::sqrt(k * m));
    double Q = (zeta <= 0.0) ? std::numeric_limits<double>::infinity() : 1.0 / (2.0 * zeta);
    return {wn, zeta, Q};
}

bool resonance_frequency(double wn, double zeta, double &wr_out) {
    const double crit = 1.0 / std::sqrt(2.0);
    if (zeta >= crit) return false;
    wr_out = wn * std::sqrt(1.0 - 2.0 * zeta * zeta);
    return true;
}

// Half-power frequencies around the resonant peak of normalized low-pass second order
bool half_power_frequencies(double wn, double zeta, double &w1, double &w2) {
    double wr;
    if (!resonance_frequency(wn, zeta, wr)) return false;
    double inside = 1.0 - zeta * zeta;
    if (inside <= 0.0) return false;

    double r1_sq = 1.0 - 2.0 * zeta * zeta - 2.0 * zeta * std::sqrt(inside);
    double r2_sq = 1.0 - 2.0 * zeta * zeta + 2.0 * zeta * std::sqrt(inside);
    if (r1_sq <= 0.0 || r2_sq <= 0.0) return false;

    w1 = wn * std::sqrt(r1_sq);
    w2 = wn * std::sqrt(r2_sq);
    return true;
}

// Force->displacement transfer function: G(s) = 1 / (m s^2 + c s + k)
std::complex<double> G_force_to_disp(std::complex<double> s, double m, double c, double k) {
    return 1.0 / (m*s*s + c*s + k);
}

int main() {
    // Example mass-spring-damper
    double m = 1.0, c = 0.4, k = 100.0;
    auto p = second_order_from_mck(m, c, k);

    std::cout << "Mass–spring–damper parameters\n";
    std::cout << "  wn   = " << p.wn   << " rad/s\n";
    std::cout << "  zeta = " << p.zeta << "\n";
    std::cout << "  Q    = " << p.Q    << "\n";

    double wr;
    if (resonance_frequency(p.wn, p.zeta, wr)) {
        std::cout << "  wr   = " << wr << " rad/s\n";
        double w1, w2;
        if (half_power_frequencies(p.wn, p.zeta, w1, w2)) {
            double bw = w2 - w1;
            std::cout << "  w1   = " << w1 << " rad/s\n";
            std::cout << "  w2   = " << w2 << " rad/s\n";
            std::cout << "  BW   = " << bw << " rad/s\n";
            std::cout << "  Q_hp = " << (wr / bw) << " (wr/BW)\n";
            std::cout << "  BW approx (2 zeta wn) = " << (2.0*p.zeta*p.wn) << " rad/s\n";
        }
    } else {
        std::cout << "  wr   = (no resonant peak; zeta >= 1/sqrt(2))\n";
    }

    // Numeric sweep to confirm peak of |G(jw)| for force->displacement model
    std::vector<double> w;
    w.reserve(2000);
    for (int i = 0; i < 2000; ++i) {
        // logspace from 1e-1 to 1e3
        double a = -1.0 + 4.0 * (double)i / (2000.0 - 1.0);
        w.push_back(std::pow(10.0, a));
    }
    double w_peak = w[0];
    double mag_peak = 0.0;
    for (double wi : w) {
        std::complex<double> s(0.0, wi);
        double mag = std::abs(G_force_to_disp(s, m, c, k));
        if (mag > mag_peak) { mag_peak = mag; w_peak = wi; }
    }
    std::cout << "\nNumerical sweep (force->displacement) peak:\n";
    std::cout << "  w_peak ~ " << w_peak << " rad/s\n";
    std::cout << "  |G(jw_peak)| ~ " << mag_peak << "\n";

    // Series RLC quick computation
    double R = 10.0, L = 50e-3, C = 10e-6;
    double w0 = 1.0 / std::sqrt(L*C);
    double Q = (w0 * L) / R;
    double BW = R / L;
    std::cout << "\nSeries RLC:\n";
    std::cout << "  w0 = " << w0 << " rad/s\n";
    std::cout << "  Q  = " << Q  << "\n";
    std::cout << "  BW = " << BW << " rad/s\n";

    return 0;
}
      

Chapter12_Lesson3.java

In Java, basic computations can be done with arrays; for complex arithmetic and numerical fitting, Apache Commons Math is a standard scientific library.

// Chapter12_Lesson3.java
// System Dynamics — Chapter 12, Lesson 3
// Resonance, Bandwidth, and Quality Factor in Mechanical and Electrical Systems
//
// Compile:
//   javac Chapter12_Lesson3.java
// Run:
//   java Chapter12_Lesson3
//
// Notes:
// - This file avoids external dependencies by using closed-form magnitude formulas.
// - For richer numerical computing (complex arithmetic, fitting), consider Apache Commons Math.

public class Chapter12_Lesson3 {

    static class SecondOrderParams {
        final double wn;   // rad/s
        final double zeta; // damping ratio
        final double Q;    // quality factor

        SecondOrderParams(double wn, double zeta, double Q) {
            this.wn = wn; this.zeta = zeta; this.Q = Q;
        }
    }

    static SecondOrderParams secondOrderFromMCK(double m, double c, double k) {
        if (m <= 0.0 || k <= 0.0) throw new IllegalArgumentException("m and k must be positive.");
        double wn = Math.sqrt(k / m);
        double zeta = c / (2.0 * Math.sqrt(k * m));
        double Q = (zeta <= 0.0) ? Double.POSITIVE_INFINITY : 1.0 / (2.0 * zeta);
        return new SecondOrderParams(wn, zeta, Q);
    }

    static Double resonanceFrequency(double wn, double zeta) {
        double crit = 1.0 / Math.sqrt(2.0);
        if (zeta >= crit) return null;
        return wn * Math.sqrt(1.0 - 2.0 * zeta * zeta);
    }

    static double[] halfPowerFrequencies(double wn, double zeta) {
        Double wr = resonanceFrequency(wn, zeta);
        if (wr == null) return null;
        double inside = 1.0 - zeta * zeta;
        if (inside <= 0.0) return null;

        double r1sq = 1.0 - 2.0*zeta*zeta - 2.0*zeta*Math.sqrt(inside);
        double r2sq = 1.0 - 2.0*zeta*zeta + 2.0*zeta*Math.sqrt(inside);
        if (r1sq <= 0.0 || r2sq <= 0.0) return null;

        return new double[] { wn*Math.sqrt(r1sq), wn*Math.sqrt(r2sq) };
    }

    // Normalized low-pass magnitude squared:
    // |G(jw)|^2 = 1 / ((1-r^2)^2 + (2 zeta r)^2), r = w/wn
    static double magNormalizedLowpass(double w, double wn, double zeta) {
        double r = w / wn;
        double denom = (1.0 - r*r)*(1.0 - r*r) + (2.0*zeta*r)*(2.0*zeta*r);
        return 1.0 / Math.sqrt(denom);
    }

    public static void main(String[] args) {
        double m = 1.0, c = 0.4, k = 100.0;
        SecondOrderParams p = secondOrderFromMCK(m, c, k);

        System.out.println("Mass–spring–damper parameters:");
        System.out.println("  wn   = " + p.wn + " rad/s");
        System.out.println("  zeta = " + p.zeta);
        System.out.println("  Q    = " + p.Q);

        Double wr = resonanceFrequency(p.wn, p.zeta);
        if (wr != null) {
            System.out.println("  wr   = " + wr + " rad/s");
            double[] hp = halfPowerFrequencies(p.wn, p.zeta);
            if (hp != null) {
                double w1 = hp[0], w2 = hp[1];
                double bw = w2 - w1;
                System.out.println("  w1   = " + w1 + " rad/s");
                System.out.println("  w2   = " + w2 + " rad/s");
                System.out.println("  BW   = " + bw + " rad/s");
                System.out.println("  Q_hp = " + (wr / bw) + " (wr/BW)");
            }
        } else {
            System.out.println("  wr   = (no resonant peak; zeta >= 1/sqrt(2))");
        }

        // Series RLC
        double R = 10.0, L = 50e-3, C = 10e-6;
        double w0 = 1.0 / Math.sqrt(L*C);
        double Q = (w0 * L) / R;
        double BW = R / L;
        System.out.println("\nSeries RLC:");
        System.out.println("  w0 = " + w0 + " rad/s");
        System.out.println("  Q  = " + Q);
        System.out.println("  BW = " + BW + " rad/s");
    }
}
      

Chapter12_Lesson3.m

MATLAB’s Control System Toolbox provides tf, bode, and related utilities. The script below also programmatically builds a small Simulink model to support frequency-sweep experiments.

% Chapter12_Lesson3.m
% System Dynamics — Chapter 12, Lesson 3
% Resonance, Bandwidth, and Quality Factor in Mechanical and Electrical Systems
%
% Toolboxes:
%   - Control System Toolbox (tf, bode)
%   - Simulink (optional section builds a small model programmatically)

clear; clc;

%% 1) Mass–spring–damper example: m xdd + c xd + k x = f
m = 1.0;     % kg
c = 0.4;     % N*s/m
k = 100.0;   % N/m

wn   = sqrt(k/m);
zeta = c/(2*sqrt(k*m));
Q    = 1/(2*zeta);

fprintf("Mass–spring–damper:\n");
fprintf("  wn   = %.6g rad/s\n", wn);
fprintf("  zeta = %.6g\n", zeta);
fprintf("  Q    = %.6g\n", Q);

% Resonant frequency exists only if zeta < 1/sqrt(2)
if zeta < 1/sqrt(2)
    wr = wn*sqrt(1 - 2*zeta^2);
    fprintf("  wr   = %.6g rad/s\n", wr);

    % Half-power frequencies relative to peak (exact)
    r1sq = 1 - 2*zeta^2 - 2*zeta*sqrt(1 - zeta^2);
    r2sq = 1 - 2*zeta^2 + 2*zeta*sqrt(1 - zeta^2);
    w1 = wn*sqrt(r1sq);
    w2 = wn*sqrt(r2sq);
    BW = w2 - w1;
    fprintf("  w1   = %.6g rad/s\n", w1);
    fprintf("  w2   = %.6g rad/s\n", w2);
    fprintf("  BW   = %.6g rad/s\n", BW);
    fprintf("  Q_hp = %.6g (wr/BW)\n", wr/BW);
    fprintf("  BW approx (2 zeta wn) = %.6g rad/s\n", 2*zeta*wn);
else
    fprintf("  wr   = (no resonant peak; zeta >= 1/sqrt(2))\n");
end

%% 2) Frequency response check using transfer functions
% Force -> displacement: G(s) = 1/(m s^2 + c s + k)
sys = tf(1, [m c k]);

w = logspace(-1, 3, 4000);  % rad/s
[mag, phase] = bode(sys, w);
mag = squeeze(mag);

% Estimate peak and half-power points numerically
[Mr, idx] = max(mag);
wr_hat = w(idx);
target = Mr / sqrt(2);

% left crossing
w1_hat = NaN;
for i = idx:-1:2
    if (mag(i)-target)*(mag(i-1)-target) <= 0
        t = (target - mag(i-1))/(mag(i) - mag(i-1));
        w1_hat = w(i-1) + t*(w(i) - w(i-1));
        break;
    end
end

% right crossing
w2_hat = NaN;
for i = idx:1:(numel(w)-1)
    if (mag(i)-target)*(mag(i+1)-target) <= 0
        t = (target - mag(i))/(mag(i+1) - mag(i));
        w2_hat = w(i) + t*(w(i+1) - w(i));
        break;
    end
end

fprintf("\nEstimated from bode samples:\n");
fprintf("  wr_hat  = %.6g rad/s\n", wr_hat);
fprintf("  Mr_hat  = %.6g\n", Mr);
fprintf("  w1_hat  = %.6g rad/s\n", w1_hat);
fprintf("  w2_hat  = %.6g rad/s\n", w2_hat);
fprintf("  BW_hat  = %.6g rad/s\n", w2_hat - w1_hat);
fprintf("  Q_hat   = %.6g\n", wr_hat/(w2_hat - w1_hat));

%% 3) Series RLC quick computation (theoretical)
R = 10.0; L = 50e-3; C = 10e-6;
w0 = 1/sqrt(L*C);
Q_rlc = w0*L/R;
BW_rlc = R/L;
fprintf("\nSeries RLC:\n");
fprintf("  w0 = %.6g rad/s\n", w0);
fprintf("  Q  = %.6g\n", Q_rlc);
fprintf("  BW = %.6g rad/s\n", BW_rlc);

%% 4) Optional: Build a small Simulink model programmatically
% The model: Sine Wave -> Transfer Fcn (1/(m s^2 + c s + k)) -> Scope
% This is a minimal setup that can be used for frequency sweeps.

build_simulink = true;
if build_simulink
    mdl = "Chapter12_Lesson3_Simulink";
    if bdIsLoaded(mdl)
        close_system(mdl, 0);
    end
    new_system(mdl); open_system(mdl);

    add_block("simulink/Sources/Sine Wave", mdl + "/Sine");
    add_block("simulink/Continuous/Transfer Fcn", mdl + "/Plant");
    add_block("simulink/Sinks/Scope", mdl + "/Scope");

    set_param(mdl + "/Plant", "Numerator", "1", "Denominator", sprintf("[%g %g %g]", m, c, k));

    add_line(mdl, "Sine/1", "Plant/1");
    add_line(mdl, "Plant/1", "Scope/1");

    set_param(mdl, "StopTime", "10");
    save_system(mdl);
    fprintf("\nCreated Simulink model: %s\n", mdl);
end
      

Chapter12_Lesson3.nb

Wolfram Mathematica provides symbolic and numeric tools for transfer functions, resonance conditions, and plots.

(* Chapter12_Lesson3.nb
   System Dynamics — Chapter 12, Lesson 3
   Resonance, Bandwidth, and Quality Factor in Mechanical and Electrical Systems
*)

Notebook[{
  Cell["Chapter 12, Lesson 3 — Resonance, Bandwidth, and Quality Factor", "Title"],
  Cell["1) Canonical second-order system (normalized low-pass)", "Section"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      G[s_] := (wn^2)/(s^2 + 2 zeta wn s + wn^2)
    ]
  ], "Input"],
  Cell["Magnitude squared and resonance condition", "Text"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      (* |G(j w)|^2 = 1 / ((1 - r^2)^2 + (2 zeta r)^2), r = w/wn *)
      M2[r_] := 1/((1 - r^2)^2 + (2 zeta r)^2);
      (* resonance when d/dr denominator = 0 => r^2 = 1 - 2 zeta^2, if zeta < 1/Sqrt[2] *)
    ]
  ], "Input"],

  Cell["2) Example: mass–spring–damper parameters", "Section"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      m = 1.0; c = 0.4; k = 100.0;
      wn = Sqrt[k/m];
      zeta = c/(2 Sqrt[k m]);
      Q = 1/(2 zeta);
      wr = If[zeta < 1/Sqrt[2], wn Sqrt[1 - 2 zeta^2], Missing["NoPeak"]];
      {wn, zeta, Q, wr}
    ]
  ], "Input"],

  Cell["Half-power frequencies relative to peak (exact)", "Text"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      r1sq = 1 - 2 zeta^2 - 2 zeta Sqrt[1 - zeta^2];
      r2sq = 1 - 2 zeta^2 + 2 zeta Sqrt[1 - zeta^2];
      w1 = wn Sqrt[r1sq];
      w2 = wn Sqrt[r2sq];
      bw = w2 - w1;
      Qhp = wr/bw;
      {w1, w2, bw, Qhp}
    ]
  ], "Input"],

  Cell["3) Plot normalized magnitude", "Section"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      Mag[w_] := 1/Sqrt[(1 - (w/wn)^2)^2 + (2 zeta (w/wn))^2];
      LogLinearPlot[Mag[w], {w, 0.1, 1000}, PlotRange -> All,
        AxesLabel -> {"w (rad/s)", "|G(j w)|"}]
    ]
  ], "Input"],

  Cell["4) Series RLC: theoretical Q and bandwidth", "Section"],
  Cell[BoxData@ToBoxes[
    HoldForm[
      R = 10.0; L = 50*10^-3; C = 10*10^-6;
      w0 = 1/Sqrt[L C];
      Qrlc = w0 L/R;
      BWrlc = R/L;
      {w0, Qrlc, BWrlc}
    ]
  ], "Input"]
}]
      

Problem 1 (Resonance condition for second order): For \( G(s) = \omega_n^2/(s^2 + 2\zeta\omega_n s + \omega_n^2) \):, prove that an interior resonant peak exists iff \( \zeta < 1/\sqrt{2} \): and derive \( \omega_r \):.

Solution: Use \( |G(j\omega)|^2 = 1/D(r) \): with \( D(r)=(1-r^2)^2+(2\zeta r)^2 \):, \( r=\omega/\omega_n \):. Minimizing \( D(r) \): gives \( dD/dr = 4r(r^2-(1-2\zeta^2)) \):. A positive critical point exists iff \( 1-2\zeta^2 > 0 \):, i.e. \( \zeta < 1/\sqrt{2} \):, and then

\[ r_r = \sqrt{1-2\zeta^2}, \quad \omega_r = \omega_n r_r = \omega_n\sqrt{1-2\zeta^2}. \]

Problem 2 (Half-power bandwidth and the approximation BW ≈ 2ζωn): For the same \( G(s) \):, derive the exact half-power frequencies relative to the resonant peak and show that for \( \zeta \ll 1 \):, \( \mathrm{BW} \approx 2\zeta\omega_n \):.

Solution: The half-power condition is \( D(r) = 2D_{\min} \):, where \( D_{\min} = 4\zeta^2(1-\zeta^2) \):. Setting \( y=r^2 \): yields a quadratic whose roots are

\[ r_{1,2}^2 = 1 - 2\zeta^2 \mp 2\zeta\sqrt{1-\zeta^2}, \quad \omega_{1,2} = \omega_n\sqrt{r_{1,2}^2}. \]

For \( \zeta \ll 1 \):, use \( \sqrt{1-\zeta^2} \approx 1 - \zeta^2/2 \): and \( \sqrt{1+\epsilon} \approx 1+\epsilon/2 \):. Then \( \omega_1 \approx \omega_n(1-\zeta) \): and \( \omega_2 \approx \omega_n(1+\zeta) \):, hence \( \mathrm{BW}=\omega_2-\omega_1 \approx 2\zeta\omega_n \):.

Problem 3 (Mechanical resonance magnitude): For \( m\ddot x + c\dot x + kx = f(t) \):, compute the resonant magnitude \( |X(j\omega_r)/F(j\omega_r)| \): in terms of \( k \): and \( \zeta \):, assuming \( \zeta < 1/\sqrt{2} \):.

Solution: Write \( X/F = (1/k)\,\omega_n^2/(s^2+2\zeta\omega_n s+\omega_n^2) \):. The normalized peak is \( 1/(2\zeta\sqrt{1-\zeta^2}) \):, so

\[ \left|\frac{X(j\omega_r)}{F(j\omega_r)}\right| = \frac{1}{k}\,\frac{1}{2\zeta\sqrt{1-\zeta^2}}. \]

Problem 4 (Exact bandwidth of a series RLC): For a series RLC excited by a sinusoidal voltage source, show that the half-power bandwidth (defined using current amplitude) is exactly \( \mathrm{BW} = R/L \): and therefore \( Q = \omega_0 L/R \):.

Solution: With \( Z(j\omega)=R + j(\omega L - 1/(\omega C)) \):, the current magnitude is \( |I| = |V|/\sqrt{R^2 + (\omega L - 1/(\omega C))^2} \):. At resonance \( \omega_0 = 1/\sqrt{LC} \):, the peak is \( |I_0|=|V|/R \):. Half-power means \( |I| = |I_0|/\sqrt{2} \):, hence \( (\omega L - 1/(\omega C))^2 = R^2 \):. Multiply by \( \omega^2 \): to obtain the two equations \( L\omega^2 - 1/C = \pm R\omega \):, i.e.

\[ L\omega^2 - R\omega - \frac{1}{C}=0 \quad \text{and} \quad L\omega^2 + R\omega - \frac{1}{C}=0. \]

The positive roots are \( \omega_2 = (R+\sqrt{R^2+4L/C})/(2L) \): and \( \omega_1 = (-R+\sqrt{R^2+4L/C})/(2L) \):. Therefore \( \mathrm{BW}=\omega_2-\omega_1=R/L \): exactly, and \( Q=\omega_0/\mathrm{BW}=\omega_0 L/R \):.

Problem 5 (Estimating ζ from Bode data): Suppose you measure a resonant peak at \( \omega_r = 50 \): rad/s and half-power points at \( \omega_1 = 48 \): rad/s and \( \omega_2 = 52 \): rad/s. Estimate \( Q \): and \( \zeta \): (light damping approximation).

Solution: The bandwidth is \( \mathrm{BW}=\omega_2-\omega_1=4 \): rad/s, hence

\[ Q \approx \frac{\omega_r}{\mathrm{BW}} = \frac{50}{4} = 12.5, \quad \zeta \approx \frac{1}{2Q} = \frac{1}{25} = 0.04. \]

If higher accuracy is required, one should use the exact half-power formulas in Section 3 to solve for \( \zeta \):.

We derived precise resonance conditions for second-order LTI systems, obtained closed-form half-power bandwidth expressions, and unified the frequency-response and energy definitions of quality factor. Mechanical mass–spring–damper and electrical RLC networks were shown to share the same second-order structure, with exact bandwidth/Q relationships in canonical RLC cases. Practical workflows and multi-language code were provided for computing and estimating \( \omega_r \):, bandwidth, and \( Q \): from theoretical models and sampled frequency-response data.