The double pendulum

Two point masses on two light rods, hanging one below the other. This is the classic example of a system that is simple to describe, coupled through inertia, and chaotic at large amplitudes. The paper derives the nonlinear and linearized equations with the Lagrangian method and simulates both in Simulink. Its potential energy misses one term, so the upper pendulum feels too little gravity. The interactive simulator reproduces the paper's model, the corrected one, normal modes and chaos.

Problem statement

A double pendulum with rod lengths $L_1$, $L_2$ and masses $m_1$, $m_2$ is shown in the figure. Determine the equations of motion of the system for arbitrary oscillations, and also for small oscillations.

The rods are rigid and massless, and the joints are frictionless.
The motion takes place in the vertical plane of the drawing.
The generalized coordinates are the absolute angles $\theta_1$ and $\theta_2$ of the two rods, both measured from the vertical.

There are no non-conservative elements, so all dynamic properties are contained in the Lagrangian. The system has two degrees of freedom.

Figure 1. Double pendulum: massless rods L₁, L₂, point masses m₁, m₂, absolute angles θ₁, θ₂ from the vertical

Energies of the system

TerminologyAs in the course: attractive energy $\mathcal{A}$ = potential energy, repulsive coenergy $\mathcal{R}^*$ = kinetic coenergy, $\mathcal{L}=\mathcal{R}^*-\mathcal{A}$.

Kinetic coenergy

Mass $m_1$ moves on a circle, so $v_1 = L_1\dot\theta_1$. The velocity of $m_2$ is the vector sum of $m_1$'s velocity $L_1\dot\theta_1$ (perpendicular to rod 1) and the relative velocity $L_2\dot\theta_2$ (perpendicular to rod 2). The angle between these two vectors is $\theta_2-\theta_1$, so the cosine rule gives

paper$$v_2^2 = (L_1\dot\theta_1)^2 + (L_2\dot\theta_2)^2 + 2L_1L_2\dot\theta_1\dot\theta_2\cos(\theta_2-\theta_1)$$

$$\mathcal{R}^* = \tfrac12(m_1+m_2)L_1^2\dot\theta_1^2 + \tfrac12 m_2L_2^2\dot\theta_2^2 + m_2L_1L_2\dot\theta_1\dot\theta_2\cos(\theta_2-\theta_1)$$

The same result follows from the Cartesian positions $x_2 = L_1\sin\theta_1 + L_2\sin\theta_2$, $y_2 = -L_1\cos\theta_1 - L_2\cos\theta_2$, differentiated and squared. This route is a useful check.

Potential energy, as written in the paper

The paper measures the height of each mass from its lowest position and writes

paper$$z_1 = L_1(1-\cos\theta_1),\qquad z_2 = L_2(1-\cos\theta_2),\qquad \mathcal{A} = m_1gz_1 + m_2gz_2$$

Missing term$z_2$ is the height of $m_2$ relative to $m_1$ only. Mass $m_2$ is carried up together with $m_1$, so its true height is $$z_2 = L_1(1-\cos\theta_1) + L_2(1-\cos\theta_2).$$ The term $m_2gL_1(1-\cos\theta_1)$ is missing, so in the paper's model the upper rod feels gravity from $m_1$ only, not from both masses. See section 7.

Lagrangian

The paper's Lagrangian (1):

(1)$$\mathcal{L} = \tfrac12 L_1^2(m_1+m_2)\dot\theta_1^2 + \tfrac12 m_2L_2^2\dot\theta_2^2 + m_2L_1L_2\dot\theta_1\dot\theta_2\cos(\theta_2-\theta_1) - m_1gL_1(1-\cos\theta_1) - m_2gL_2(1-\cos\theta_2)$$

The corrected Lagrangian has $(m_1+m_2)$ in the first gravity term:

(1′)$$\mathcal{L} = \tfrac12 (m_1+m_2)L_1^2\dot\theta_1^2 + \tfrac12 m_2L_2^2\dot\theta_2^2 + m_2L_1L_2\dot\theta_1\dot\theta_2\cos(\theta_2-\theta_1) - (m_1+m_2)gL_1(1-\cos\theta_1) - m_2gL_2(1-\cos\theta_2)$$

Equations of motion (arbitrary oscillations)

With no non-conservative forces, Lagrange's equation for each coordinate is

(2)$$\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot\theta_i}\right) - \frac{\partial\mathcal{L}}{\partial\theta_i} = 0,\qquad i = 1,2$$

Coordinate $\theta_1$

$$\frac{\partial\mathcal{L}}{\partial\dot\theta_1} = L_1^2(m_1+m_2)\dot\theta_1 + m_2L_1L_2\dot\theta_2\cos(\theta_2-\theta_1),\qquad \frac{\partial\mathcal{L}}{\partial\theta_1} = m_2L_1L_2\dot\theta_1\dot\theta_2\sin(\theta_2-\theta_1) - m_1gL_1\sin\theta_1$$

Differentiating in time, the two $\dot\theta_1\dot\theta_2\sin$ terms cancel. The paper shows this cancellation explicitly. The result is

(*)$$(m_1+m_2)L_1^2\ddot\theta_1 + m_2L_1L_2\ddot\theta_2\cos(\theta_2-\theta_1) - m_2L_1L_2\dot\theta_2^2\sin(\theta_2-\theta_1) + m_1gL_1\sin\theta_1 = 0$$

Coordinate $\theta_2$

$$\frac{\partial\mathcal{L}}{\partial\dot\theta_2} = m_2L_2^2\dot\theta_2 + m_2L_1L_2\dot\theta_1\cos(\theta_2-\theta_1),\qquad \frac{\partial\mathcal{L}}{\partial\theta_2} = -m_2L_1L_2\dot\theta_1\dot\theta_2\sin(\theta_2-\theta_1) - m_2gL_2\sin\theta_2$$

(**)$$m_2L_2^2\ddot\theta_2 + m_2L_1L_2\ddot\theta_1\cos(\theta_2-\theta_1) + m_2L_1L_2\dot\theta_1^2\sin(\theta_2-\theta_1) + m_2gL_2\sin\theta_2 = 0$$

In the paper, $\partial\mathcal{L}/\partial\theta_2$ is written with a "+" in front of the velocity term. Since $\partial\cos(\theta_2-\theta_1)/\partial\theta_2 = -\sin(\theta_2-\theta_1)$, the sign should be "−". A second sign slip in the next line cancels this one, so the final (**) is correct.

Corrected equations in matrix form

With the missing gravity term restored, only the last term of (*) changes, from $m_1$ to $m_1+m_2$. Written as $\mathbf M(\boldsymbol\theta)\ddot{\boldsymbol\theta} = \mathbf f$, the equations are

(3)$$\begin{bmatrix}(m_1+m_2)L_1^2 & m_2L_1L_2\cos\Delta\\ m_2L_1L_2\cos\Delta & m_2L_2^2\end{bmatrix}\begin{bmatrix}\ddot\theta_1\\ \ddot\theta_2\end{bmatrix} = \begin{bmatrix} m_2L_1L_2\dot\theta_2^2\sin\Delta - (m_1+m_2)gL_1\sin\theta_1\\ -m_2L_1L_2\dot\theta_1^2\sin\Delta - m_2gL_2\sin\theta_2\end{bmatrix},\qquad \Delta = \theta_2-\theta_1$$

The mass matrix is symmetric, as it must be, and its determinant $m_2L_1^2L_2^2(m_1 + m_2\sin^2\Delta)$ is always positive. It can therefore always be inverted, and this is how the simulator below computes $\ddot\theta_1$ and $\ddot\theta_2$.

Small oscillations

The paper sets $\sin\theta\approx\theta$, $\cos\theta\approx1$, $\cos(\theta_2-\theta_1)\approx1$ and $\sin(\theta_2-\theta_1)\approx0$, and obtains

(*-1)$$(m_1+m_2)L_1^2\ddot\theta_1 + m_2L_1L_2\ddot\theta_2 - m_2L_1L_2\dot\theta_2^2 + m_1gL_1\theta_1 = 0$$

(**-1)$$m_2L_2^2\ddot\theta_2 + m_2L_1L_2\ddot\theta_1 + m_2L_1L_2\dot\theta_1^2(\theta_2-\theta_1) + m_2gL_2\theta_2 = 0$$

These are not linear. In (*-1) the factor $\sin(\theta_2-\theta_1)$ was dropped but the term $\dot\theta_2^2$ kept, which contradicts the stated approximation. A consistent linearization keeps only terms of first order in the angles and angular velocities. Both velocity-squared terms are of second or third order and vanish:

(4)$$\begin{aligned}(m_1+m_2)L_1^2\ddot\theta_1 + m_2L_1L_2\ddot\theta_2 + (m_1+m_2)gL_1\theta_1 &= 0\\ m_2L_1L_2\ddot\theta_1 + m_2L_2^2\ddot\theta_2 + m_2gL_2\theta_2 &= 0\end{aligned}$$

or $\mathbf M\ddot{\boldsymbol\theta} + \mathbf K\boldsymbol\theta = \mathbf 0$, with a constant symmetric $\mathbf M$ and the diagonal $\mathbf K = \operatorname{diag}\big((m_1+m_2)gL_1,\; m_2gL_2\big)$. The pendulums are coupled only through inertia (the off-diagonal $m_2L_1L_2$), not through stiffness. This contrasts with pendulums joined by a spring, as in assignment 10.

Simulink models

Nonlinear model (equations (*) and (**)): two chains of two integrators, one for $\theta_1$ and one for $\theta_2$. A subtract block forms $\theta_2-\theta_1$, which feeds sin and cos blocks. Product blocks form $\cos(\theta_2-\theta_1)\,\ddot\theta_j$ and $\sin(\theta_2-\theta_1)\,\dot\theta_j^2$, and gains scale each term. The gravity terms come from sin blocks on $\theta_1$ and $\theta_2$. Parameters: $m_1=0.5$, $m_2=1$, $L_1=L_2=1$, $g=9.81$.
Small-oscillation model (equations (*-1) and (**-1)): the same structure without trigonometric blocks, with two product blocks for the velocity-squared terms. Parameters: $m_1=m_2=1$, $L_1=L_2=1$, $g=9.81$.

Algebraic loopEach equation contains the other angle's acceleration ($\ddot\theta_2$ in (*), $\ddot\theta_1$ in (**)). Wiring an integrator input back into the other chain creates an algebraic loop. Simulink can solve it iteratively at each step, but this is slower and can fail. The cleaner approach is to solve the 2×2 system (3) for both accelerations first, e.g. in a MATLAB Function block. The gain values are not shown in the paper, so the diagram cannot be fully checked.

Simulation results from the paper

Plot	Model and parameters	Initial state (read from the plot)
Plot 1	Nonlinear, $m_1=0.5$, $m_2=1$	θ₁ = θ₂ ≈ 0.07 rad (4°), at rest
Plot 2	Small oscillations, $m_1=m_2=1$	θ₁ = θ₂ ≈ 1 rad (57°), at rest

Both plots show the two angles over time, in yellow and magenta, without a legend or time labels. The curves are visibly polygonal, because Simulink's default output step is coarse. The paper concludes that "for arbitrary and for small oscillations the motions differ very little from each other".

The comparison is not validThe two runs use different masses ($m_1 = 0.5$ vs 1) and different amplitudes (4° vs 57°). Two different systems in two different regimes cannot show whether linearization is accurate. With equal parameters and equal initial states (presets "Fair 4°" and "Fair 57°" in the simulator), the linear and nonlinear models agree to within 0.001 rad at 4°. At 57° they differ by more than 1.5 rad within 10 s.

We were not able to reproduce the plots closely with either the paper's or the corrected equations. The time scale is unknown and the Simulink gains are not given. The overall shape, two frequencies mixed with a slow and a faster component, is typical of every double pendulum.

Critical review

The potential energy of $m_2$ is incomplete

The term $m_2gL_1(1-\cos\theta_1)$ is missing, so in the paper's equations the gravity torque on the upper rod is $m_1gL_1\sin\theta_1$ instead of $(m_1+m_2)gL_1\sin\theta_1$. One check exposes this immediately: if $m_1\to0$, the paper's upper rod has no restoring torque at all, even though a heavy $m_2$ hangs below it. The error changes every frequency:

Parameters ($L_1=L_2=1$)	Paper's model ω₁, ω₂ [rad/s]	Correct model ω₁, ω₂ [rad/s]
$m_1=0.5,\ m_2=1$ (plot 1)	1.62, 6.05	2.32, 7.31
$m_1=m_2=1$ (plot 2)	1.94, 5.07	2.40, 5.79

The slow mode is up to 30 % too slow. Even at 4°, the angles from the paper's model differ from the correct ones by up to 0.17 rad within 10 s, which is more than the amplitude itself. The paper's model also does not conserve the true mechanical energy (see the energy chart in the simulator).

The "small-oscillation" equations are not linear

(*-1) keeps $-m_2L_1L_2\dot\theta_2^2$ after setting $\sin(\theta_2-\theta_1)=0$. This term is quadratic, never negative, and pushes $\theta_1$ in one direction only. (**-1) keeps $\dot\theta_1^2(\theta_2-\theta_1)$, which is of third order. Both should vanish, leaving the linear system (4). Only a truly linear system has normal modes and a closed-form solution.

"Small" oscillations simulated at 57°

Plot 2 starts from about 1 rad. That is far outside the range where $\sin\theta\approx\theta$ holds: at 57°, $\sin\theta$ is already 16 % smaller than $\theta$. The linear model is meant for a few degrees.

Smaller points

Equation (2) mixes indices $i$ and $j$ (it should be $i$ in both terms). The equations of motion are written without "= 0". The plots lack a legend, axis labels with units and a time scale. A finer output step (or the Refine option) would give smooth curves. The paper says nothing about the most interesting property of the double pendulum: at large amplitudes its motion is chaotic (section 8).

Extension: normal modes and chaos

Normal modes

For the linear system (4), $\det(\mathbf K-\omega^2\mathbf M) = 0$. For equal lengths $L_1=L_2=L$ and the mass ratio $\mu = m_2/m_1$, this gives

(5)$$\omega_{1,2}^2 = \frac gL\Big[(1+\mu) \mp \sqrt{\mu(1+\mu)}\Big]$$

For equal masses ($\mu=1$), $\omega^2 = (2\mp\sqrt2)\,g/L$. The mode shapes are $\theta_2/\theta_1 = +\sqrt2$ (slow, both rods swing to the same side) and $\theta_2/\theta_1 = -\sqrt2$ (fast, rods move in opposite directions). Starting exactly in a mode shape gives a pure single-frequency oscillation. The presets "Mode 1" and "Mode 2" demonstrate this. The paper's initial state $\theta_1=\theta_2$ is a mixture of both modes, which is why its plots show two frequencies.

Chaos

For small amplitudes the motion is regular and quasi-periodic. With enough energy (e.g. both rods released from 90° or more), the lower rod can flip over, and the motion becomes chaotic. Two runs that start 0.001° apart separate completely after a few seconds:

θ₁(0) = θ₂(0)	Max difference in θ₂ after 30 s (initial difference 0.001°)
10°	0.00003 rad (regular)
60°	0.0001 rad (regular)
90°	≈ 14 rad (chaotic, separation after ≈ 15 s)
120°	≈ 48 rad (chaotic, separation after ≈ 8 s)

Equal masses and lengths are assumed. In the chaotic regime no simulation, Simulink or otherwise, can predict the long-term motion, because tiny rounding errors grow in the same way. Energy is still conserved, which makes it a good check of the numerical solution.

Interactive simulator

The simulator solves the equations numerically (4th-order Runge–Kutta, 1 ms step), with the 2×2 mass matrix inverted at every step. The paper's models are included, so its results can be compared with the corrected ones.

Cases

Model

Compare with

δ for perturbed run [°]

System parameters

m₁ [kg]

m₂ [kg]

L₁ [m]

L₂ [m]

g [m/s²]

θ₁(0) [°]

θ₂(0) [°]

θ̇₁(0) [rad/s]

θ̇₂(0) [rad/s]

Duration [s]

Angle fields accept expressions, e.g. 180/pi for 1 rad.

Angles θ₁(t) and θ₂(t)

Difference between the two runs

True mechanical energy

Evaluated with the correct energy expression. Constant only for the correct nonlinear model.

Path of the lower mass m₂

Animation (selected model)

Results table (every 0.05 s)

Review questions

How can you detect the missing gravity term without re-deriving the energy?

Use a limiting case. If $m_1\to0$, the system is a single mass $m_2$ on a jointed rod, and gravity must still pull the upper rod back to the vertical. The paper's equation (*) has no gravity term left in that case. Another check: in static equilibrium, the upper joint carries the weight of both masses.

Why are the pendulums coupled even though there is no spring between them?

The motion of the upper mass moves the pivot of the lower pendulum. In the equations this shows up as the off-diagonal term $m_2L_1L_2\cos(\theta_2-\theta_1)$ of the mass matrix: inertial coupling. The stiffness matrix of the linear model is diagonal.

Why must the term $\dot\theta_2^2$ be dropped in the small-oscillation equations?

Linearization keeps terms of first order in the small quantities $\theta_i$ and $\dot\theta_i$. $\dot\theta_2^2\sin(\theta_2-\theta_1)$ is of third order. Dropping only the sine and keeping $\dot\theta_2^2$ produces a second-order term that does not belong to either the linear or the full model.

How would you start the pendulum so that it oscillates at a single frequency?

Release it from rest in a mode shape, an eigenvector of $\mathbf M^{-1}\mathbf K$. For equal masses and lengths: $\theta_2 = \sqrt2\,\theta_1$ (slow mode) or $\theta_2 = -\sqrt2\,\theta_1$ (fast mode), with small amplitudes.

What does "chaotic" mean for a deterministic system like this one?

The equations have no randomness, but two nearby initial states separate exponentially fast. Any error in the initial state or in the numerical integration is amplified, so the long-term motion cannot be predicted, although the total energy stays constant.

How would you avoid the algebraic loop in Simulink?

Solve equation (3) for $\ddot\theta_1$ and $\ddot\theta_2$ explicitly (Cramer's rule, or a MATLAB Function block computing $\mathbf M^{-1}\mathbf f$), and feed each result directly into its own integrator chain.

Explanation of Method