Forced Oscillations and Resonance

Learning Objectives

By the end of this section you will be able to:

1. Write the driven oscillator equation $m\ddot{x} + c\dot{x} + kx = F_0\cos(\omega t)$ and identify which of its terms produce the transient and which produce the steady state.
2. Derive the steady-state amplitude $A(\omega)$ and phase lag $\varphi(\omega)$ using the complex-exponential (impedance) method.
3. Explain resonance physically — what happens to a lightly damped system when the driving frequency approaches the natural frequency.
4. Compute the resonant frequency, the maximum amplitude, and the quality factor $Q$.
5. Recognize beats in the undamped limit and the unbounded resonance solution $x_p(t) = \tfrac{F_0}{2m\omega_0}\,t\sin(\omega_0 t)$.
6. Simulate a driven oscillator in Python (RK4) and sweep frequencies in PyTorch with one complex line.

The Driving Question

“A child on a swing weighs 20 kg. You weigh 80 kg. Why can the child push the swing higher than you can?”

Anyone who has watched a playground knows the answer: the child pushes at the right rhythm. Every nudge arrives just as the swing reaches its forward edge, adding energy stride for stride. An adult applying ten times the force at the wrong rhythm gets nowhere — the push sometimes accelerates, sometimes brakes, averaging out.

That intuition is resonance. A small, repeated, perfectly-timed input causes a huge response. The mathematics of this section is just the precise language to say it. By the end you'll know exactly which rhythm, exactly how big the response gets, and exactly why it grows.

The Equation of a Driven Oscillator

Take the mass–spring–damper from Section 22.06 and add an external force $F(t)$ pushing the mass. Newton's second law says

The most important and most studied case is a sinusoidal force at a single frequency $\omega$:

Transient + Steady State

Because the equation is linear, its general solution splits into two pieces:

Finding the Steady-State Solution

We need a particular solution. The classical guess is $x_p = A\cos(\omega t - \varphi)$, but the clean way to derive it uses complex exponentials. Replace the real forcing with $F_0 e^{i\omega t}$ (we'll take the real part at the end) and guess

Differentiating twice multiplies by $i\omega$, so $\dot{x}_p = i\omega X e^{i\omega t}$ and $\ddot{x}_p = -\omega^2 X e^{i\omega t}$. Plugging in:

Cancel the common factor and solve for $X$:

The denominator is called the (mechanical) impedance $Z(\omega)$, in direct analogy with the electrical impedance $Z(\omega) = R + i\omega L + \tfrac{1}{i\omega C}$ of an RLC circuit. One complex number per frequency contains everything: take its modulus to get the amplitude, its argument to get the phase.

Amplitude and Phase Lag — The Core Formulas

Write $X = |X|e^{-i\varphi}$ in polar form (negative sign so $\varphi$ is the phase the mass lags the force by). Taking the real part of $X e^{i\omega t}$ recovers the real steady state:

where the amplitude and phase are obtained from $|Z(\omega)|$ and $\arg Z(\omega)$:

These two formulas are the heart of every vibration textbook in existence. Memorize them — or better, learn to re-derive them in 30 seconds with the impedance trick above.

Resonance: Why a Small Push Can Become Huge

Look at the amplitude formula $A(\omega) = F_0/\sqrt{(k-m\omega^2)^2 + (c\omega)^2}$ and pretend you can pick $\omega$. You want $A$ to be as large as possible, so you want the denominator to be as small as possible.

The first term $(k - m\omega^2)^2$ can be made exactly zero by choosing $\omega = \sqrt{k/m} = \omega_0$. Then only the damping term $(c\omega_0)^2$ is left, and

If $c \to 0$, this blows up. That is the mathematical story behind every “a bridge fell down because the wind oscillated at the bridge's natural frequency” anecdote.

The quality factor $Q$ packages the “sharpness” of the peak into one number:

At resonance $A(\omega_0) = F_0/(c\omega_0) = (F_0/k)\cdot Q$: the static deflection $F_0/k$ is amplified by Q. A car tire ($Q \approx 5$), a guitar string ($Q \approx 1000$), a crystal oscillator ($Q \approx 10^6$) and a LIGO mirror ($Q \approx 10^7$) span twelve orders of magnitude on the same single equation.

Interactive: The Universal Resonance Curve

Below is the dimensionless resonance curve $A(r) = 1/\sqrt{(1-r^2)^2 + (2\zeta r)^2}$ with $r = \omega/\omega_0$. Slide $\zeta$ from 0.05 (sharp peak) toward 1 (peak disappears) and watch the bandwidth widen as the peak shrinks. The phase plot beneath it always crosses $\pi/2$ at $r = 1$ — that is the universal signature of resonance.

Interactive: Live Driven Oscillator

Now watch the whole story unfold in time. Start with the default values ($m=1, c=0.4, k=4, F_0=2, \omega=2$) — the system is being driven exactly at its natural frequency $\omega_0 = \sqrt{k/m} = 2$. The amber curve starts at zero and builds toward the cyan steady-state envelope at $A = 2.5$. Then try sliding $\omega$ off the resonance — the amplitude collapses almost immediately.

Worked Example (Step-by-Step)

Let's do every step by hand for the system in the simulator: $m = 1$, $c = 0.4$, $k = 4$, $F_0 = 2$, $\omega = 2$, $x(0) = 0$, $x'(0) = 0$.

The Undamped Case: Beats and Pure Resonance

Set $c = 0$. The amplitude formula becomes $A(\omega) = F_0/|k - m\omega^2|$, which is finite for every $\omega \neq \omega_0$ but blows up at $\omega = \omega_0$. What actually happens?

Off-resonance: beats

When $\omega$ is close to but not equal to $\omega_0$, with initial conditions zero:

Using the sum-to-product identity $\cos A - \cos B = -2\sin\!\bigl(\tfrac{A+B}{2}\bigr)\sin\!\bigl(\tfrac{A-B}{2}\bigr)$,

The fast factor at frequency $(\omega_0+\omega)/2$ looks like a normal oscillation, but it's amplitude-modulated by the slow factor at frequency $(\omega_0-\omega)/2$ — that's a beat. You hear beats when you tune a guitar string near a reference pitch: the volume swells and fades at the difference frequency.

On resonance: secular (linear-in-t) growth

The above formula has a 0/0 problem at $\omega = \omega_0$. Take the limit carefully (L'Hôpital, or use the method of undetermined coefficients with $x_p = A t \sin(\omega_0 t)$) and you find

The amplitude $\tfrac{F_0}{2m\omega_0}\,t$ grows linearly forever. In an ideal undamped system, a perfectly tuned force pumps unlimited energy into the oscillator. Real systems are saved by either nonzero damping (the system never quite gets there) or by nonlinearity (the spring stops being Hookean at large displacements).

Computation in Python

Two ingredients: the algebraic $A(\omega), \varphi(\omega)$ formulas, and a plain RK4 integrator that solves the full ODE. We check that the integrator's late-time output matches the algebra.

1import numpy as np
2import matplotlib.pyplot as plt
3
4def steady_state(m, c, k, F0, omega):
5    """
6    Closed-form amplitude and phase lag of
7        m x'' + c x' + k x = F0 cos(omega t).
8    Returns (A, phi) so that  x_ss(t) = A cos(omega t - phi).
9    """
10    denom_re = k - m * omega**2
11    denom_im = c * omega
12    Z = np.sqrt(denom_re**2 + denom_im**2)        # |impedance|
13    A   = F0 / Z
14    phi = np.arctan2(denom_im, denom_re)          # in [0, pi]
15    return A, phi
16
17def simulate(m, c, k, F0, omega, x0, v0, t):
18    """RK4 integration of the FULL driven equation."""
19    A_sys = np.array([[0.0,    1.0],
20                      [-k/m, -c/m]])
21    def rhs(tt, Y):
22        force = np.array([0.0, F0 * np.cos(omega * tt) / m])
23        return A_sys @ Y + force
24
25    Y = np.array([x0, v0], dtype=float)
26    out = np.zeros((len(t), 2))
27    out[0] = Y
28    for i in range(1, len(t)):
29        h  = t[i] - t[i - 1]
30        k1 = rhs(t[i - 1],         Y)
31        k2 = rhs(t[i - 1] + h/2,   Y + h * k1 / 2)
32        k3 = rhs(t[i - 1] + h/2,   Y + h * k2 / 2)
33        k4 = rhs(t[i - 1] + h,     Y + h * k3)
34        Y  = Y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
35        out[i] = Y
36    return out[:, 0]
37
38# ---- the worked example ----
39m, c, k = 1.0, 0.4, 4.0
40F0, omega = 2.0, 2.0              # drive AT the natural frequency
41t = np.linspace(0, 30, 1500)
42
43A, phi = steady_state(m, c, k, F0, omega)
44print(f"steady-state amplitude  A   = {A:.4f}")     # 2.5000
45print(f"phase lag               phi = {np.degrees(phi):.2f} deg")  # 90 deg
46
47x_num = simulate(m, c, k, F0, omega, x0=0.0, v0=0.0, t=t)
48x_ss  = A * np.cos(omega * t - phi)
49
50print(f"|x_num(end) - x_ss(end)| = {abs(x_num[-1] - x_ss[-1]):.2e}")  # ~0

PyTorch: Complex-Impedance Sweep

The complex form $X(\omega) = F_0 / Z(\omega)$ is one line of vectorized math. PyTorch's native complex tensors let us sweep hundreds of frequencies at once, then pull out amplitude, phase, and quality factor with a few more lines. This is exactly how vibration engineers compute “frequency-response functions” in production tools.

1import torch
2
3def resonance_sweep(m, c, k, F0, omegas):
4    """
5    Sweep driving frequencies and return complex steady-state amplitude X(omega).
6
7        m x'' + c x' + k x = F0 e^{i omega t}    (complex form)
8        guess  x = X e^{i omega t}
9        =>  (-m omega^2 + i c omega + k) X = F0
10        =>  X(omega) = F0 / ( k - m omega^2  +  i c omega )
11    """
12    re = k - m * omegas**2
13    im = c * omegas
14    Z  = torch.complex(re, im)            # impedance, one complex per frequency
15    return F0 / Z                          # complex amplitude X(omega)
16
17# physical setup
18m, c, k = 1.0, 0.4, 4.0
19F0      = 2.0
20omega_0 = (k / m) ** 0.5                  # natural frequency
21
22omegas = torch.linspace(0.05, 4.0, 400)
23X       = resonance_sweep(m, c, k, F0, omegas)
24
25A   = X.abs()                             # amplitude  |X|
26phi = -X.angle()                          # phase lag  (mass lags force by phi)
27
28# locate the resonance peak
29peak_idx = torch.argmax(A)
30print(f"natural freq        omega_0    = {omega_0:.4f}")
31print(f"peak amplitude      A_max      = {A[peak_idx].item():.4f}")
32print(f"peak frequency      omega_peak = {omegas[peak_idx].item():.4f}")
33
34# Q factor and bandwidth (half-power points)
35A_half = A[peak_idx] / (2.0 ** 0.5)
36mask   = A >= A_half
37bw     = omegas[mask].max() - omegas[mask].min()
38print(f"half-power bandwidth Delta_omega = {bw.item():.4f}")
39print(f"quality factor       Q          = {(omega_0 / bw).item():.3f}")

Where Resonance Shows Up

Every row uses the same $m\ddot{x} + c\dot{x} + kx = F_0\cos(\omega t)$ equation. The only thing changing is the units of $m, c, k, x$. That is the power of having one canonical form: solve it once, recognize it everywhere.

Common Pitfalls

Summary