Coupled oscillators and normal modes

Hang two identical pendulums side by side and join them with a weak spring. Pull one aside and let go. It swings, but its swings slowly die away while the other pendulum, at first motionless, starts swinging more and more. After a while the first pendulum is still and the second is swinging fully; then the energy flows back. This slow exchange looks complicated, but it is the sum of two very simple motions, the normal modes, in which both pendulums swing at a single common frequency. This lesson shows how to find normal modes for any set of coupled oscillators, why every small motion is a superposition of them, and how a long chain of coupled masses turns into a continuous string with a wave equation. Normal modes are the bridge from oscillations to waves, and the same mathematics describes molecular vibrations, the tones of a drum and the sound modes of a crystal.

Two coupled pendulums

Let each pendulum have mass \(m\) and length \(\ell\), so on its own it has \(\omega_0^2 = g/\ell\). Join the bobs with a spring of stiffness \(k\), unstretched when both hang vertically, and use the horizontal displacements \(x_1, x_2\) as coordinates. For small swings the gravitational restoring force on each bob is \(-m\omega_0^2 x\), and the spring is stretched by \(x_2 - x_1\), pulling bob 1 toward bob 2 with force \(k(x_2 - x_1)\). So

The equations are coupled: each acceleration depends on both coordinates. Add and subtract them. With \(q_1 = (x_1 + x_2)/2\) and \(q_2 = (x_1 - x_2)/2\), and writing \(\kappa = k/m\),

In these normal coordinates the equations decouple into two independent harmonic oscillators. The first mode, \(q_2 = 0\), has \(x_1 = x_2\): the pendulums swing together, the spring never stretches, and the frequency is \(\omega_1 = \omega_0\). The second mode, \(q_1 = 0\), has \(x_1 = -x_2\): the pendulums swing in opposite directions, the spring is stretched by twice the displacement of each, and the frequency is \(\omega_2 = \sqrt{\omega_0^2 + 2\kappa}\), higher because the spring adds stiffness. Every motion of the system is a combination \(x_1 = q_1 + q_2\), \(x_2 = q_1 - q_2\) of the two modes with arbitrary amplitudes and phases.

Energy exchange as beating

Start with \(x_1 = A\), \(x_2 = 0\) and both at rest. Then \(q_1 = q_2 = A/2\) initially, so each mode is excited with amplitude \(A/2\) and

For weak coupling, \(\omega_2 - \omega_1\) is small, so each pendulum swings at nearly \(\omega_0\) inside a slowly varying envelope. The envelope of \(x_1\) falls to zero just as that of \(x_2\) reaches its maximum, after a time \(\pi/(\omega_2 - \omega_1)\). For \(\kappa \ll \omega_0^2\), \(\omega_2 - \omega_1 \approx \kappa/\omega_0\), so the transfer time is about \(\pi\omega_0/\kappa\): the weaker the spring, the slower the exchange. The same beating appears in sound when two tuning forks of slightly different pitch sound together, and in quantum mechanics as the oscillation of an electron between two coupled wells.

The general method

For any number of coupled oscillators near equilibrium, the kinetic and potential energies are quadratic forms, \(T = \tfrac12\dot{\mathbf x}^{\mathsf T}M\dot{\mathbf x}\) and \(U = \tfrac12\mathbf x^{\mathsf T}K\mathbf x\), with a symmetric positive-definite mass matrix \(M\) and a symmetric stiffness matrix \(K\). The equations of motion are \(M\ddot{\mathbf x} = -K\mathbf x\). Look for a normal mode, a motion in which every coordinate oscillates at the same frequency with fixed relative amplitudes, \(\mathbf x = \mathbf a\cos(\omega t + \phi)\). Substituting gives the generalized eigenvalue problem

The determinant is a polynomial of degree \(n\) in \(\omega^2\), so there are \(n\) values \(\omega_j^2\), each with a mode shape \(\mathbf a_j\). Because \(M\) and \(K\) are symmetric and \(M\) is positive definite, the eigenvalues are real, and if \(K\) is positive definite they are positive. The modes are orthogonal in the sense \(\mathbf a_i^{\mathsf T}M\mathbf a_j = 0\) for \(i \ne j\). To prove it, write \(K\mathbf a_i = \omega_i^2 M\mathbf a_i\) and \(K\mathbf a_j = \omega_j^2 M\mathbf a_j\). Multiply the first on the left by \(\mathbf a_j^{\mathsf T}\) and the second by \(\mathbf a_i^{\mathsf T}\), and subtract. By symmetry the left sides are equal, so \((\omega_i^2 - \omega_j^2)\mathbf a_j^{\mathsf T}M\mathbf a_i = 0\), and distinct frequencies force the product to vanish. Orthogonality is what lets any initial condition be decomposed uniquely into modes, just as the Fourier series of lesson 7 decomposes any shape of a string.

A chain of N masses

Generalize to \(N\) equal masses joined by \(N + 1\) springs between fixed walls, with spacing \(a\). Mass \(j\) obeys \(m\ddot x_j = k(x_{j+1} - 2x_j + x_{j-1})\), with \(x_0 = x_{N+1} = 0\). Try a mode in which the amplitudes vary sinusoidally along the chain, \(x_j = \sin(j\theta)\cos\omega t\). Using \(\sin((j+1)\theta) + \sin((j-1)\theta) = 2\sin(j\theta)\cos\theta\), the equation becomes \(-m\omega^2 = k(2\cos\theta - 2)\), so

The wall condition \(x_{N+1} = 0\) requires \(\sin((N+1)\theta) = 0\), so \(\theta = n\pi/(N+1)\) for \(n = 1, \ldots, N\). These are the \(N\) normal modes, with frequencies \(\omega_n = 2\sqrt{k/m}\,\sin(n\pi/2(N+1))\). For \(N = 3\) this gives \(2\sin(22.5^\circ) = 0.765\), \(2\sin(45^\circ) = 1.414\) and \(2\sin(67.5^\circ) = 1.848\), matching the worked example.

Write \(\theta = \kappa_w a\), so that \(x_j = \sin(\kappa_w \cdot ja)\) is a sinusoid in the position \(ja\) with wavenumber \(\kappa_w\). The frequency formula \(\omega = 2\sqrt{k/m}\,|\sin(\kappa_w a/2)|\) is a dispersion relation: it says how frequency depends on wavenumber. For long wavelengths, \(\kappa_w a \ll 1\), the sine is approximately its argument and \(\omega \approx \sqrt{k/m}\,a\,\kappa_w\), proportional to \(\kappa_w\), which is the behaviour of a nondispersive wave with speed \(v = a\sqrt{k/m}\). At short wavelengths the chain is dispersive, and there is a maximum frequency \(2\sqrt{k/m}\) above which no wave propagates. A crystal behaves like this: its atoms form a lattice, sound is the long-wavelength limit of its vibrations, and the maximum frequency is a few terahertz.

The continuum limit

Now let the number of masses grow while the total length \(L = (N + 1)a\) stays fixed. Each mass \(m = \mu a\) for a linear density \(\mu\), and each spring under tension behaves, for transverse motion, like \(k = T/a\) with tension \(T\). Then \(v = a\sqrt{k/m} = \sqrt{T/\mu}\), independent of \(a\). The difference equation becomes

because the bracket is the standard finite-difference approximation to the second derivative. This is the wave equation for a string, derived directly in lesson 5. The lowest modes of the chain become the modes of the string, \(y_n = \sin(n\pi x/L)\cos\omega_n t\) with \(\omega_n = n\pi v/L\), a harmonic series. Only the highest modes, whose wavelengths are comparable to the spacing \(a\), remember the discreteness. The program checks how quickly the chain frequencies approach the string frequencies.

Normal modes in nature

A molecule with \(n\) atoms has \(3n\) coordinates and so \(3n\) modes, of which six for a nonlinear molecule (five for a linear one) are rigid translations and rotations with zero frequency. The remaining vibrational modes have frequencies in the infrared, and each absorbs light at its own frequency, which is how infrared spectroscopy identifies molecules. Carbon dioxide, a linear molecule with three atoms, has four vibrational modes: a symmetric stretch near \(1340\ \mathrm{cm^{-1}}\) in the units spectroscopists use, an antisymmetric stretch near \(2350\ \mathrm{cm^{-1}}\), and two degenerate bending modes near \(667\ \mathrm{cm^{-1}}\). The bending modes absorb strongly in the thermal infrared that Earth radiates, which is the physical basis of the greenhouse effect of carbon dioxide. Engineering structures have modes too, and a structural engineer's first job with a new tower or bridge deck is to find its lowest few natural frequencies and keep them away from those of wind, footsteps and machinery.

Computing normal modes

The program solves the generalized eigenvalue problem with a small Jacobi rotation routine written from scratch, checks the three-mass result, confirms the \(N\)-mass formula, and compares the lowest chain frequencies with the continuous string as \(N\) grows.

import math

def jacobi_eigen(A, tol=1e-12):
    """Eigenvalues and eigenvectors of a symmetric matrix by Jacobi rotations."""
    n = len(A); A = [row[:] for row in A]
    V = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
    while True:
        p, q, big = 0, 1, 0.0
        for i in range(n):
            for j in range(i + 1, n):
                if abs(A[i][j]) > big: p, q, big = i, j, abs(A[i][j])
        if big < tol: break
        th = 0.5 * math.atan2(2 * A[p][q], A[q][q] - A[p][p])
        c, s = math.cos(th), math.sin(th)
        for k in range(n):
            akp, akq = A[k][p], A[k][q]
            A[k][p], A[k][q] = c*akp - s*akq, s*akp + c*akq
        for k in range(n):
            apk, aqk = A[p][k], A[q][k]
            A[p][k], A[q][k] = c*apk - s*aqk, s*apk + c*aqk
        for k in range(n):
            vkp, vkq = V[k][p], V[k][q]
            V[k][p], V[k][q] = c*vkp - s*vkq, s*vkp + c*vkq
    vals = [A[i][i] for i in range(n)]
    order = sorted(range(n), key=lambda i: vals[i])
    return [vals[i] for i in order], [[V[k][i] for k in range(n)] for i in order]

def chain(N):            # stiffness matrix for N unit masses, unit springs, fixed ends
    return [[2.0 if i == j else (-1.0 if abs(i - j) == 1 else 0.0) for j in range(N)] for i in range(N)]

vals, vecs = jacobi_eigen(chain(3))
print("three masses: omega / omega_s =", [round(math.sqrt(v), 4) for v in vals])
for v in vecs:
    s = v[0] if abs(v[0]) > 1e-9 else 1.0
    print("  mode shape", [round(c / s, 4) + 0.0 for c in v])

N = 8
vals, _ = jacobi_eigen(chain(N))
formula = [2 * math.sin(n * math.pi / (2 * (N + 1))) for n in range(1, N + 1)]
print(f"\nN = {N}: largest difference from the formula = {max(abs(math.sqrt(a) - b) for a, b in zip(vals, formula)):.2e}")

# continuum limit: string of length 1, wave speed 1, so omega_n = n pi
print("\nlowest three chain frequencies (string limit pi, 2pi, 3pi):")
for N in (3, 10, 30, 100):
    a = 1.0 / (N + 1)
    w = [2 / a * math.sin(n * math.pi * a / 2) for n in (1, 2, 3)]
    print(f"N = {N:3d}: " + ", ".join(f"{x:.4f}" for x in w))

# coupled pendulums of the worked example
w0, kap = math.sqrt(9.81), 0.40
w2 = math.sqrt(w0*w0 + 2*kap)
print(f"\ncoupled pendulums: w1 = {w0:.4f}, w2 = {w2:.4f}, transfer time = {math.pi/(w2 - w0):.2f} s")

three masses: omega / omega_s = [0.7654, 1.4142, 1.8478]
  mode shape [1.0, 1.4142, 1.0]
  mode shape [1.0, 0.0, -1.0]
  mode shape [1.0, -1.4142, 1.0]

N = 8: largest difference from the formula = 1.11e-15

lowest three chain frequencies (string limit pi, 2pi, 3pi):
N =   3: 3.0615, 5.6569, 7.3910
N =  10: 3.1309, 6.1981, 9.1391
N =  30: 3.1402, 6.2724, 9.3885
N = 100: 3.1415, 6.2822, 9.4214

coupled pendulums: w1 = 3.1321, w2 = 3.2573, transfer time = 25.09 s

The Jacobi routine reproduces the three frequencies and mode shapes of the worked example and agrees with the \(N\)-mass formula to rounding error. In the continuum check the chain frequencies, scaled to a string of unit length and unit wave speed, approach \(\pi\), \(2\pi\) and \(3\pi\) as \(N\) grows. Already at \(N = 30\) the fundamental is within a small fraction of a part in a hundred, while the higher modes converge more slowly because their wavelengths are closer to the spacing.