Simple harmonic motion and energy

Lesson 1 showed that any system near a smooth stable equilibrium behaves like a mass on a spring. This lesson solves that problem completely. We derive the general solution of the equation of motion, fit it to initial conditions, follow the energy as it passes back and forth between kinetic and potential forms, draw the motion in phase space, and meet the same equation in a pendulum swinging about a pivot, a liquid sloshing in a U-tube, a torsion balance and an electrical circuit. We finish by integrating the equation numerically with a method that respects energy conservation, the velocity Verlet algorithm, which is the workhorse of molecular dynamics and orbital simulation.

The equation and its solutions

A mass \(m\) on a spring of stiffness \(k\), sliding without friction, obeys Newton's second law \(m\ddot x = -kx\), or

This is a linear second-order differential equation with constant coefficients. Try \(x = e^{st}\): substituting gives \(s^2 + \omega_0^2 = 0\), so \(s = \pm i\omega_0\). The general solution is a combination \(x = C_1 e^{i\omega_0 t} + C_2 e^{-i\omega_0 t}\), and for \(x\) to be real we need \(C_2 = \bar C_1\). Writing \(C_1 = \tfrac12 A e^{i\phi}\) gives the familiar real form

with \(a = A\cos\phi\) and \(b = -A\sin\phi\). Two constants are exactly what a second-order equation needs: they are fixed by the initial position \(x_0\) and velocity \(v_0\). Setting \(t = 0\) in \(x\) and in \(\dot x = -a\omega_0\sin\omega_0 t + b\omega_0\cos\omega_0 t\) gives \(a = x_0\) and \(b = v_0/\omega_0\), so

Uniqueness matters here. The existence and uniqueness theorem for linear differential equations guarantees that the solution with given \(x_0\) and \(v_0\) is the only one, so once we have found it by any method we are done. The velocity \(\dot x = -A\omega_0\sin(\omega_0 t + \phi)\) leads the displacement by a quarter cycle, and the acceleration \(\ddot x = -\omega_0^2 x\) is exactly opposite the displacement, with magnitude largest at the turning points where the velocity is zero.

There is a picture behind the formula. A point moving at constant angular speed \(\omega_0\) around a circle of radius \(A\) has an \(x\) coordinate \(A\cos(\omega_0 t + \phi)\). Simple harmonic motion is the shadow of uniform circular motion, and the phase \(\omega_0 t + \phi\) is the angle around the circle. This is the phasor of lesson 1 made literal: the complex number \(Ae^{i(\omega_0 t + \phi)}\) rotates in the complex plane and its real part is the motion.

Energy

Multiply the equation of motion by \(\dot x\): \(m\dot x\ddot x + kx\dot x = 0\), which is \(\frac{d}{dt}\big(\tfrac12 m\dot x^2 + \tfrac12 kx^2\big) = 0\). The total energy

is constant. Substituting the solution, the kinetic energy is \(\tfrac12 kA^2\sin^2(\omega_0 t + \phi)\) and the potential energy is \(\tfrac12 kA^2\cos^2(\omega_0 t + \phi)\). Each oscillates at twice the frequency of the motion, between zero and \(E\), and they are exactly out of step. Averaged over a cycle each is \(E/2\), because the average of \(\sin^2\) or \(\cos^2\) over a cycle is one half. This equal sharing is special to the quadratic potential; for a potential \(U \propto |x|^n\) the virial theorem gives \(\langle K \rangle = \tfrac n2 \langle U \rangle\), which reduces to equality for \(n = 2\).

Energy also gives the speed at any position without solving for the time dependence: \(v = \pm\omega_0\sqrt{A^2 - x^2}\). The cart in the worked example moves at \(15\sqrt{0.0640^2 - 0.050^2} = 0.60\ \mathrm{m/s}\) at \(x = 5.0\ \mathrm{cm}\), as it should, and at the maximum speed \(\omega_0 A = 0.96\ \mathrm{m/s}\) as it crosses equilibrium.

Phase space

The state of the oscillator at any instant is the pair \((x, \dot x)\), a point in the phase plane. As time passes the point moves around a closed curve. From the energy equation, \(\tfrac12 m\dot x^2 + \tfrac12 kx^2 = E\) is an ellipse with semi-axes \(A\) along \(x\) and \(\omega_0 A\) along \(\dot x\). Plotting \(\dot x/\omega_0\) instead of \(\dot x\) turns the ellipses into circles traversed clockwise at angular speed \(\omega_0\). Larger energies give larger curves, and curves for different energies never cross, because the state at one instant determines the whole future.

For the pendulum the phase portrait shows what the harmonic approximation misses. Small orbits are nearly circles. Larger ones stretch out horizontally, and at energy \(2mg\ell\) they meet the separatrix, a curve through the unstable upright position \(\theta = \pm\pi\) that divides swinging from rotating motion. Above it the pendulum goes over the top and the orbits are wavy lines rather than loops. The phase portrait is the best single picture of what a nonlinear oscillator can do.

The same equation in other systems

A physical pendulum is any rigid body swinging about a pivot. If its mass is \(m\), its moment of inertia about the pivot is \(I\), and its centre of mass is a distance \(d\) below the pivot, the gravitational torque for a small angle is \(-mgd\theta\), and \(I\ddot\theta = -mgd\,\theta\) gives \(\omega_0 = \sqrt{mgd/I}\). A uniform rod of length \(\ell\) pivoted at one end has \(I = \tfrac13 m\ell^2\) and \(d = \ell/2\), so \(\omega_0 = \sqrt{3g/2\ell}\), faster than a simple pendulum of the same length because its mass is spread closer to the pivot.

A torsion pendulum, a disc hanging from a wire that twists, obeys \(I\ddot\theta = -\kappa\theta\) with the torsion constant \(\kappa\) of the wire. Henry Cavendish used one in 1798 to measure the gravitational attraction between lead spheres. Liquid in a U-tube with total length \(\ell\) of liquid column oscillates with \(\omega_0 = \sqrt{2g/\ell}\), because displacing the surfaces by \(x\) creates a height difference \(2x\). A floating object with cross-sectional area \(A\) bobs with \(\omega_0 = \sqrt{\rho g A / m}\). And an inductor \(L\) in series with a capacitor \(C\) obeys \(L\ddot q + q/C = 0\), so the charge oscillates at \(\omega_0 = 1/\sqrt{LC}\), with inductance playing the role of mass and inverse capacitance the role of stiffness.

Beyond the small-angle limit

For a pendulum the next term in the potential, \(-\theta^4/24\) relative to \(\theta^2/2\), softens the restoring torque at large angles, so the pendulum spends longer near its turning points and the period increases. A perturbation calculation, or the series expansion of the exact elliptic-integral result, gives

At \(\theta_0 = 0.5\) rad, about \(29^\circ\), the correction is \(1 + 0.0156 + 0.0002\), which matches the exact ratio 1.0159. Huygens found the cure for this in 1673: a bob that moves along a cycloid instead of a circle has a period exactly independent of amplitude, a property called tautochronism. In practice clockmakers kept the amplitude small and constant instead.

Integrating the equation numerically

Most oscillators of practical interest have no closed-form solution, so we integrate numerically. The obvious method, Euler's, updates \(x \gets x + v\,\Delta t\) and \(v \gets v + a\,\Delta t\), but for an oscillator it adds energy at every step, so the orbit spirals outward. The velocity Verlet method does much better:

It is time-reversible and symplectic, which means it preserves phase-space area exactly. As a result its energy error stays bounded and oscillates rather than drifting, even over millions of steps. The program integrates a harmonic oscillator with both methods, compares the energy after 100 periods, and then integrates a pendulum at \(60^\circ\) amplitude to measure its period from the zero crossings.

import math

def euler(x, v, acc, dt, n):
    for _ in range(n):
        x, v = x + v * dt, v + acc(x) * dt
    return x, v

def verlet(x, v, acc, dt, n):
    a = acc(x)
    for _ in range(n):
        x = x + v * dt + 0.5 * a * dt * dt
        a_new = acc(x)
        v = v + 0.5 * (a + a_new) * dt
        a = a_new
    return x, v

w0 = 2 * math.pi            # period 1 s
acc = lambda x: -w0 * w0 * x
energy = lambda x, v: 0.5 * v * v + 0.5 * w0 * w0 * x * x
E0 = energy(1.0, 0.0)
dt, steps = 0.01, 100 * 100   # 100 periods
for name, method in (("Euler", euler), ("Verlet", verlet)):
    x, v = method(1.0, 0.0, acc, dt, steps)
    print(f"{name:6s}: x = {x:+.5f}, energy ratio after 100 periods = {energy(x, v)/E0:.6f}")

# pendulum at 60 degrees: measure the period from upward zero crossings of theta
g, L = 9.81, 1.0
acc_p = lambda th: -(g / L) * math.sin(th)
th, om, t, dt = math.radians(60), 0.0, 0.0, 1e-4
a = acc_p(th); crossings = []
while len(crossings) < 3:
    th_new = th + om * dt + 0.5 * a * dt * dt
    a_new = acc_p(th_new)
    om = om + 0.5 * (a + a_new) * dt
    if th < 0 <= th_new:                       # crossing zero going up
        crossings.append(t + dt * (-th) / (th_new - th))
    th, a, t = th_new, a_new, t + dt
T_num = crossings[2] - crossings[1]
T0 = 2 * math.pi * math.sqrt(L / g)
th0 = math.radians(60)
print(f"pendulum at 60 deg: T = {T_num:.5f} s, T/T0 = {T_num/T0:.5f}")
print(f"series 1 + th0^2/16 + 11 th0^4/3072 = {1 + th0**2/16 + 11*th0**4/3072:.5f}")

Euler : x = +244002576.66156, energy ratio after 100 periods = 129271963281139616.000000
Verlet: x = +0.99466, energy ratio after 100 periods = 0.999989
pendulum at 60 deg: T = 2.15287 s, T/T0 = 1.07318
series 1 + th0^2/16 + 11 th0^4/3072 = 1.07285

Euler's method multiplies the energy by \(1 + \omega_0^2\Delta t^2\) each step, so after 10,000 steps the energy has grown by a factor of about \(10^{17}\); the cart would have flown off the track. Verlet keeps the energy within about one part in a hundred thousand. For the pendulum at \(60^\circ\) the measured period ratio 1.0732 agrees with the exact value 1.07318 from lesson 1, and the three-term series gives 1.0729, already close.