Central forces and rigid bodies revisited with Lagrangians

PHYS-244 spent three lessons on central forces and one on rigid bodies, deriving everything from vector equations. This lesson revisits both with the Lagrangian method and shows how much of that work becomes routine: the reduction of the two-body problem, the effective potential, and the conserved angular momentum all drop out of the Lagrangian once the right coordinates are chosen. We then treat two systems that are hard with vectors and natural with Lagrangians: the spherical pendulum, whose orbit precesses, and the heavy symmetric top, whose axis precesses and nods.

The two-body problem in one stroke

Two particles of masses \(m_1\), \(m_2\) interact through a potential \(V(|\mathbf r_1 - \mathbf r_2|)\). Change to the center-of-mass position \(\mathbf R = (m_1\mathbf r_1 + m_2\mathbf r_2)/M\) and the relative position \(\mathbf r = \mathbf r_1 - \mathbf r_2\), with \(M = m_1 + m_2\) and reduced mass \(\mu = m_1m_2/M\). The kinetic energy separates:

\(\mathbf R\) is cyclic in all three components, so the total momentum \(M\dot{\mathbf R}\) is conserved and the center of mass moves uniformly. What remains is a single body of mass \(\mu\) in a fixed central potential. Rotational symmetry means \(\mathbf L = \mu\mathbf r\times\dot{\mathbf r}\) is conserved, so the motion lies in the plane perpendicular to \(\mathbf L\). In that plane, with polar coordinates,

The angle \(\phi\) is cyclic, so \(\ell = \mu r^2\dot\phi\) is conserved. The radial equation is \(\mu\ddot r - \mu r\dot\phi^2 + V'(r) = 0\), and substituting \(\dot\phi = \ell/\mu r^2\) gives

The radial motion is one-dimensional motion in the effective potential \(U_{\text{eff}} = V + \ell^2/2\mu r^2\), whose second term is the centrifugal barrier. That is the whole of PHYS-244 lesson 7 in a paragraph. The pieces were the same; the Lagrangian simply organized them.

A caution about substitution. It is tempting to substitute \(\dot\phi = \ell/\mu r^2\) into the Lagrangian itself, getting \(\tfrac12\mu\dot r^2 + \ell^2/2\mu r^2 - V\), and to read off an effective potential from it. That gives the centrifugal term with the wrong sign. Conserved momenta must be eliminated in the equations of motion, or through the Routhian (a function that replaces \(L\) by \(L - \ell\dot\phi\)), not by substituting into \(L\).

The spherical pendulum

A bob of mass \(m\) hangs from a string of length \(\ell\) and can swing in any direction. Use spherical angles: \(\theta\) from the downward vertical and azimuth \(\phi\). The bob moves with speed components \(\ell\dot\theta\) and \(\ell\sin\theta\,\dot\phi\), and its height below the pivot is \(\ell\cos\theta\), so

\(\phi\) is cyclic: \(p_\phi = m\ell^2\sin^2\theta\,\dot\phi\) is conserved, the angular momentum about the vertical. The \(\theta\) equation is

A special solution has constant \(\theta = \theta_0\): the conical pendulum, with \(\dot\phi^2 = g/(\ell\cos\theta_0)\). A slightly disturbed conical pendulum wobbles in \(\theta\) about \(\theta_0\) at a frequency \(\omega_\theta\) found from the curvature of the effective potential \(U_{\text{eff}}(\theta) = p_\phi^2/(2m\ell^2\sin^2\theta) - mg\ell\cos\theta\). If \(\omega_\theta\) were exactly twice the spin rate, the orbit seen from above would be a fixed ellipse, closest to the vertical twice per revolution. It is not, so the ellipse precesses.

Euler angles

A rigid body has six degrees of freedom: three for the position of a reference point and three for its orientation. Orientation is described by the three Euler angles \((\phi, \theta, \psi)\). Starting from a fixed frame, rotate by \(\phi\) about the vertical \(z\) axis, then by \(\theta\) about the new \(x\) axis (the line of nodes), then by \(\psi\) about the body's own symmetry axis. For a top, \(\theta\) is the tilt of the axis from the vertical, \(\phi\) is the azimuth of the axis (its precession angle), and \(\psi\) is the spin about the axis.

The angular velocity components along the body's principal axes are combinations of \(\dot\phi\), \(\dot\theta\) and \(\dot\psi\). For a symmetric body, with moments of inertia \(I_1 = I_2\) about axes perpendicular to the symmetry axis and \(I_3\) about it, the kinetic energy simplifies to

The first term is the energy of the axis swinging, exactly like the spherical pendulum; the second is the energy of spin about the axis, where \(\omega_3 = \dot\psi + \cos\theta\,\dot\phi\) is the total angular velocity about the symmetry axis.

The heavy symmetric top

A top spins with its tip fixed on a table and its center of mass a distance \(l\) up its axis. Its potential energy is \(Mgl\cos\theta\), so

with \(I_1\) measured about the tip. This Lagrangian, first solved by Lagrange, has two cyclic coordinates and no explicit time, so it has three conserved quantities:

\(p_\psi\) is the spin angular momentum, \(p_\phi\) the angular momentum about the vertical. Solving the second for \(\dot\phi\) and substituting into the energy leaves a single equation for \(\theta\), motion in an effective potential, just like the spherical pendulum. The axis therefore nods between two tilt angles while it precesses, a motion called nutation.

For a fast top, released from rest with its axis tilted, there is a simple approximate result. The gravitational torque \(Mgl\sin\theta\) acts horizontally, perpendicular to the spin angular momentum \(I_3\omega_3\), and swings it around the vertical at the rate

Superimposed on this slow precession is a fast, small nutation at a frequency close to \(I_3\omega_3/I_1\), whose amplitude shrinks as the spin grows. The program checks both.

import math

def rk4(f, y, dt):
    k1 = f(y); k2 = f([a + dt / 2 * b for a, b in zip(y, k1)])
    k3 = f([a + dt / 2 * b for a, b in zip(y, k2)]); k4 = f([a + dt * b for a, b in zip(y, k3)])
    return [a + dt / 6 * (p + 2 * q + 2 * r + s) for a, p, q, r, s in zip(y, k1, k2, k3, k4)]
g, l = 9.81, 1.0

# 1. Spherical pendulum: p_phi = m l^2 sin^2(th) phidot is conserved; th obeys
#    thddot = sin th cos th phidot^2 - (g/l) sin th
th0 = math.radians(30)
w_conic = math.sqrt(g / (l * math.cos(th0)))            # conical pendulum spin rate
h = math.sin(th0) ** 2 * w_conic                         # p_phi / (m l^2)
# small oscillations about the cone: U_eff = h^2/(2 sin^2) - (g/l) cos, frequency^2 = U_eff''
def Ueff(t): return h * h / (2 * math.sin(t) ** 2) - g / l * math.cos(t)
d = 1e-4; w_th = math.sqrt((Ueff(th0 + d) - 2 * Ueff(th0) + Ueff(th0 - d)) / d ** 2)
print(f"conical pendulum at 30 deg: spin {w_conic:.4f} rad/s; small wobble frequency {w_th:.4f} rad/s")
print(f"predicted azimuth advance between successive lowest points: {360 * w_conic / w_th:.2f} deg")
def f(y):                                                # state (th, phi, thdot)
    th, ph, thd = y; phd = h / math.sin(th) ** 2
    return [thd, phd, math.sin(th) * math.cos(th) * phd * phd - g / l * math.sin(th)]
y = [th0 + 0.02, 0.0, 0.0]; dt = 1e-4; lows = []; prev = y[2]
while len(lows) < 3:
    y = rk4(f, y, dt)
    if prev < 0 <= y[2]: lows.append(math.degrees(y[1]))     # th is lowest when thdot crosses 0 upward
    prev = y[2]
print(f"numerical azimuth advance: {lows[2] - lows[1]:.2f} deg")

# 2. Heavy symmetric top: steady precession rate for a fast top is about M g l / (I3 omega3)
I1, I3, Mgl, w3 = 8e-5, 4e-5, 0.1 * 9.81 * 0.04, 400.0
p_psi = I3 * w3
th = math.radians(30); p_phi = p_psi * math.cos(th)           # released with phidot = 0
def top(y):
    th, ph, thd = y
    phd = (p_phi - p_psi * math.cos(th)) / (I1 * math.sin(th) ** 2)
    thdd = (I1 * phd * phd * math.sin(th) * math.cos(th) - p_psi * phd * math.sin(th) + Mgl * math.sin(th)) / I1
    return [thd, phd, thdd]
y = [th, 0.0, 0.0]; dt = 1e-5; T = 2.0; thmax = th
for _ in range(int(T / dt)):
    y = rk4(top, y, dt); thmax = max(thmax, y[0])
print(f"fast-top estimate M g l / (I3 omega3) = {Mgl / (I3 * w3):.3f} rad/s; average precession rate {y[1] / T:.3f} rad/s")
print(f"nutation: the tilt from vertical swings between 30.000 and {math.degrees(thmax):.3f} deg; nutation frequency about I3 omega3 / I1 = {I3 * w3 / I1:.0f} rad/s")

conical pendulum at 30 deg: spin 3.3657 rad/s; small wobble frequency 6.0675 rad/s
predicted azimuth advance between successive lowest points: 199.69 deg
numerical azimuth advance: 199.71 deg
fast-top estimate M g l / (I3 omega3) = 2.453 rad/s; average precession rate 2.473 rad/s
nutation: the tilt from vertical swings between 30.000 and 30.726 deg; nutation frequency about I3 omega3 / I1 = 200 rad/s

A 1 m conical pendulum at 30 degrees turns at 3.366 rad/s, and small wobbles about that cone have frequency 6.068 rad/s, slightly less than twice the spin rate. The azimuth advances \(360^\circ\times3.366/6.068 = 199.7^\circ\) between successive closest approaches to the vertical, and the full integration gives 199.71 degrees. For the top, with \(I_3 = 4\times10^{-5}\ \mathrm{kg\,m^2}\), \(I_1 = 8\times10^{-5}\ \mathrm{kg\,m^2}\), \(Mgl = 0.0392\ \mathrm{N\,m}\) and spin 400 rad/s, the formula predicts precession at 2.45 rad/s; the average measured over 2 s is 2.47 rad/s. Released from rest at 30 degrees, the axis dips only to 30.73 degrees, a nutation of less than a degree at a frequency near 200 rad/s, too fast and too small to see with the naked eye. That is why a fast top seems to precess smoothly without falling.

When a top falls

As friction slows the spin, \(I_3\omega_3\) falls and the precession speeds up. Below a critical spin the top can no longer stay upright: the condition for a top standing vertically (a sleeping top) to be stable is \(I_3^2\omega_3^2 > 4I_1Mgl\). For the toy top above, the critical spin is \(\sqrt{4I_1Mgl}/I_3 = 89\ \mathrm{rad/s}\), about 850 revolutions per minute. Above that, a vertical top sleeps quietly; below it, the slightest disturbance grows and the top begins to wobble and fall.