Double pendulum

The setup

A rod of length \ell_1 hangs from a fixed pivot with a bob of mass m_1 at its end; a second rod of length \ell_2 hangs from that bob and carries m_2. The rods are rigid and massless and the bobs are points, so two angles from the vertical, \theta_1 and \theta_2, say where everything is. The lab fixes \ell_1 = 1 m and m_1 = 1 and lets you set the ratios m_2/m_1 and \ell_2/\ell_1; gravity is the standard value g = 9.806\,65 m/s², which is exact by definition [4].

The Lagrangian

With Newton's laws you'd have to track the tension in both rods. Lagrange's method lets you skip them. Write the kinetic energy T and potential energy V in terms of the two angles, and the equations of motion fall out of L = T - V. The bob positions are

x_1 = \ell_1\sin\theta_1, \quad y_1 = -\ell_1\cos\theta_1, \qquad x_2 = x_1 + \ell_2\sin\theta_2, \quad y_2 = y_1 - \ell_2\cos\theta_2,

and differentiating and adding up gives

L = \tfrac12(m_1 + m_2)\ell_1^2\dot\theta_1^2 + \tfrac12 m_2\ell_2^2\dot\theta_2^2 + m_2\ell_1\ell_2\dot\theta_1\dot\theta_2\cos(\theta_1 - \theta_2) + (m_1 + m_2)g\ell_1\cos\theta_1 + m_2 g\ell_2\cos\theta_2 .

The cross term is what ties the two pendulums together, because the lower bob's speed depends on how fast the upper one is moving. The Euler-Lagrange equation \frac{d}{dt}\frac{\partial L}{\partial\dot\theta_i} = \frac{\partial L}{\partial\theta_i} for each angle gives

(m_1 + m_2)\ell_1\ddot\theta_1 + m_2\ell_2\ddot\theta_2\cos(\theta_1 - \theta_2) + m_2\ell_2\dot\theta_2^2\sin(\theta_1 - \theta_2) + (m_1 + m_2)g\sin\theta_1 = 0,

m_2\ell_2\ddot\theta_2 + m_2\ell_1\ddot\theta_1\cos(\theta_1 - \theta_2) - m_2\ell_1\dot\theta_1^2\sin(\theta_1 - \theta_2) + m_2 g\sin\theta_2 = 0 .

That's two linear equations in the two angular accelerations. The lab solves them for \ddot\theta_1 and \ddot\theta_2, rewrites the system as four first-order equations for (\theta_1, \omega_1, \theta_2, \omega_2) and steps them forward. Its tests check that the solved accelerations make both equations above vanish. Nobody has a closed-form solution, and for most starts there isn't one to find.

Since L doesn't depend on time explicitly, the energy E = T + V is conserved. That's the lab's main check on its own numerics.

Small swings: normal modes

For small angles, \sin\theta \approx \theta, \cos\theta \approx 1 and the velocity-squared terms are negligible, so the equations become linear. With equal masses and lengths they read, in matrix form,

\begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix}\ddot\theta_1 \\ \ddot\theta_2\end{pmatrix} = -\frac{g}{\ell}\begin{pmatrix} 2 & 0 \\ 0 & 1\end{pmatrix}\begin{pmatrix}\theta_1 \\ \theta_2\end{pmatrix} .

Look for motions where both angles oscillate at one frequency, \theta = \mathbf a\cos\omega t. That is a generalised eigenvalue problem, and its two solutions are the normal modes:

• the slow mode, \theta_2 = \sqrt2\,\theta_1, with the bobs swinging together at \omega^2 = (2 - \sqrt2)\,g/\ell, a period of 2.62 s for \ell = 1 m;
• the fast mode, \theta_2 = -\sqrt2\,\theta_1, with the bobs swinging against each other at \omega^2 = (2 + \sqrt2)\,g/\ell, a period of 1.09 s.

Every small motion is a mix of the two. Nudge only the top bob and you start both, so the motion beats as energy sloshes back and forth. The ratio of the two periods is irrational, so the pattern never exactly repeats. The lab's tests run each mode for one period and check it comes back where it started.

Chaos

Pull the bobs up higher and the small-angle picture breaks down. From about 120 degrees the lower bob flips over the top at moments you can't predict, and the motion never settles into a pattern. What makes this chaos isn't how complicated it looks. It's sensitive dependence on initial conditions: two starts that differ by a tiny amount d_0 drift apart exponentially,

d(t) \approx d_0\,e^{\lambda t},

until the separation is as large as the motion itself. \lambda, the largest Lyapunov exponent, sets a horizon you can't see past. An uncertainty of 10^{-9} rad grows to order one after \ln(10^9)/\lambda \approx 21/\lambda seconds. Measure a thousand times better and you buy only 7/\lambda seconds more.

The lab runs a twin, started off by d_0 in \theta_2, and plots the distance between the two in phase space, (\theta_1, \omega_1, \theta_2, \omega_2), on a log scale. Exponential growth shows up as a straight rising line. The readout \ln(d/d_0)/t estimates \lambda while the line stays straight. Once the twins are far apart, the estimate means nothing. In the high release preset it holds near 1.2 per second. Shinbrot, Grebogi, Wisdom and Yorke measured the growth of uncertainties on a real double pendulum and found \lambda = 7.5 \pm 1.5\ \mathrm{s}^{-1}, against 7.9 \pm 0.4\ \mathrm{s}^{-1} from simulating their idealised model [1]. Their exponent is larger than the lab's because rates scale like \sqrt{g/\ell} and the lab's rods are a metre long. Levien and Tan built a version with optical encoders on the axles that quantifies the same sensitivity in an undergraduate lab [2].

Swing gently and the chaos goes away. At low energy the motion stays on smooth tori in phase space, and nearby starts drift apart only slowly. Try the low swing preset: after a minute the twins are still within 10^{-8} of each other.

Integrators and energy

The lab steps forward with the classical fourth order Runge-Kutta method. Its error per step falls like \Delta t^5, so at the default 2 ms step the energy holds to about a part in a million over many seconds. But RK4 isn't symplectic. It doesn't preserve the phase-space structure of a Hamiltonian system, so its energy error piles up instead of wobbling within a bound [3]. At a 50 ms step you can watch the energy drain away, almost half the energy scale gone in a minute. The orbit lab pits RK4 against symplectic methods on Kepler's problem, where the difference is starker.

Chaos amplifies integration error exactly as it amplifies a change in the start. Past the predictability horizon, no integrator follows the true trajectory, however small its step. What a good integrator keeps is the statistics, the energy and the shape of the motion.

Things to try in the lab

• Slow normal mode and Fast normal mode: time ten swings of each against the stated periods.
• Both modes at once: watch energy slosh between the bobs as the two modes drift in and out of step.
• Low swing, regular and High release, chaos: compare them on the log plot, one nearly flat line and one straight climb.
• Flung from a high start: the lower bob flips over the top again and again, never in the same rhythm twice.
• Step too large: the same release at a 50 ms step, with |\Delta E| climbing as RK4 drains the energy.
• Drag either bob to a new start. It starts from rest when you let go, so where you leave the bobs sets the energy.
• Make the lower bob light (mass ratio 0.1). The upper bob barely notices it and swings almost like a simple pendulum, while the lower one gets thrown around.
• Change the twin's offset from 10^{-9} to 10^{-6} and see the separation reach order one about 7/\lambda seconds sooner.

Reading the lab's instruments

Energy per kg of m_1 is T + V with the zero of potential energy at the pivot. |\Delta E| / energy scale is the drift from the start divided by (m_1 + m_2)g\ell_1 + m_2 g\ell_2, the energy needed to lift both bobs from hanging straight down to the pivot's height, a fixed scale that never passes through zero. Distance to the twin is the Euclidean distance between the two states in (\theta_1, \omega_1, \theta_2, \omega_2), with angle differences folded into one turn.

References

[1] T. Shinbrot, C. Grebogi, J. Wisdom and J. A. Yorke, "Chaos in a double pendulum," American Journal of Physics 60, 491-499 (1992). DOI: 10.1119/1.16860. The exponents are quoted from the abstract.

[2] R. B. Levien and S. M. Tan, "Double pendulum: An experiment in chaos," American Journal of Physics 61, 1038-1044 (1993). DOI: 10.1119/1.17335

[3] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration, 2nd ed. (Springer, 2006). DOI: 10.1007/3-540-30666-8

[4] NIST, CODATA Internationally Recommended Values of the Fundamental Physical Constants: standard acceleration of gravity, 9.806\,65 m s⁻², exact.