Three pendulums coupled by springs

Three identical pendulums hang side by side, and two springs connect them at mid-length. We derive the equations of motion in two ways: by the angular-momentum (torque) theorem and by Lagrange's equations. We then linearize them for small oscillations and simulate both models. Normal modes explain the striking energy exchange between the pendulums, and an interactive simulator shows where the small-angle model stops being accurate.

Problem statement

For the system in the figure, the springs are unstressed when the pendulums are vertical. Determine the equations of motion, and also consider the case of small oscillations.

Figure 1. Three pendulums of length 2a with point masses m (massless rods), pivots O₁, O₂, O₃. Springs k connect the rods at distance a from the pivots (points A, B, C). Angles α, β, γ from the vertical.

Each pendulum is a massless rod of length $2a$ with a point mass $m$ at the end. Its moment of inertia about the pivot is therefore

(12)$$I = m(2a)^2 = 4ma^2$$

First derivation: the angular-momentum theorem

The paper first works directly with vectors. The torque of a force and the angular momentum about the pivot are

(1)–(5)$$\vec M = \vec r\times\vec F,\qquad \vec L = \vec r\times m\vec v = I\vec\omega,\qquad \frac{d\vec L}{dt} = \vec r\times\vec F = \vec M\quad\Rightarrow\quad I\dot{\vec\omega} = \sum\vec M$$

Pendulum 1

Two forces act on it: gravity at the mass, and the spring force at point A. Using $\vec x\times\vec y = \vec z$:

(9)$$\vec M_1 = \overrightarrow{O_1C_1}\times m\vec g = -2amg\sin\alpha\,\vec z$$

The paper models the spring force as $k$ times the change of the vector between the attachment points, $\vec F_1 = k(\overrightarrow{AB} - \overrightarrow{O_1O_2})$. After expanding the cross products:

(10)$$\vec M_2 = \overrightarrow{O_1A}\times\vec F_1 = -a^2k\sin(\alpha-\beta)\,\vec z$$

(13)$$4ma^2\ddot\alpha + 2amg\sin\alpha + a^2k\sin(\alpha-\beta) = 0\quad\Longrightarrow\quad \ddot\alpha + \frac{g}{2a}\sin\alpha + \frac{k}{4m}\sin(\alpha-\beta) = 0$$

Pendulum 2 and pendulum 3

The middle pendulum feels both springs. The same procedure gives

(17)$$\ddot\beta + \frac{g}{2a}\sin\beta + \frac{k}{4m}\Big(\sin(\beta-\alpha) + \sin(\beta-\gamma)\Big) = 0$$

(20)$$\ddot\gamma + \frac{g}{2a}\sin\gamma + \frac{k}{4m}\sin(\gamma-\beta) = 0$$

These are three coupled nonlinear second-order equations. They are symmetric under exchanging pendulums 1 and 3.

Small oscillations

The paper justifies the approximation $\sin\theta\approx\theta$ in two ways. Geometrically, the tangent to $\sin\theta$ at $\theta=0$ has slope $\cos0 = 1$, so $g(\theta)=\theta$. Analytically, the Taylor series $\sin\theta = \theta - \theta^3/3! + \dots$ and $\cos\theta = 1 - \theta^2/2! + \dots$ can be truncated after the first term. Both give the linear system

(24)–(26)$$\begin{aligned}\ddot\alpha + \Big(\frac{g}{2a} + \frac{k}{4m}\Big)\alpha - \frac{k}{4m}\beta &= 0\\ \ddot\beta + \Big(\frac{g}{2a} + \frac{k}{2m}\Big)\beta - \frac{k}{4m}\alpha - \frac{k}{4m}\gamma &= 0\\ \ddot\gamma + \Big(\frac{g}{2a} + \frac{k}{4m}\Big)\gamma - \frac{k}{4m}\beta &= 0\end{aligned}$$

Second derivation: Lagrange's equations

The paper then notes that the problem can be solved "shorter and easier" with Lagrange's equations.

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

Each mass moves on a circle of radius $2a$, so $v_i = 2a\omega_i$:

(30)$$\mathcal{R}^* = \tfrac12 m\,4a^2(\dot\alpha^2+\dot\beta^2+\dot\gamma^2) = 2ma^2(\dot\alpha^2+\dot\beta^2+\dot\gamma^2)$$

Each mass rises by $2a(1-\cos\theta) = 4a\sin^2\tfrac\theta2$. The paper takes the spring elongations as the horizontal separations of the attachment points, $a(\sin\beta-\sin\alpha)$ and $a(\sin\gamma-\sin\beta)$:

$$\mathcal{A} = 4mga\Big(\sin^2\tfrac\alpha2 + \sin^2\tfrac\beta2 + \sin^2\tfrac\gamma2\Big) + \tfrac12k\Big[a^2(\sin\gamma-\sin\beta)^2 + a^2(\sin\beta-\sin\alpha)^2\Big]$$

For small angles:

(31)$$\mathcal{A} = mga(\alpha^2+\beta^2+\gamma^2) + \tfrac12 ka^2\Big[(\gamma-\beta)^2 + (\beta-\alpha)^2\Big]$$

(33)$$\mathcal{L} = 2ma^2(\dot\alpha^2+\dot\beta^2+\dot\gamma^2) - mga(\alpha^2+\beta^2+\gamma^2) - \tfrac12 ka^2\Big[(\gamma-\beta)^2 + (\beta-\alpha)^2\Big]$$

The Lagrange equations, for example for $\alpha$:

$$\frac{\partial\mathcal{L}}{\partial\dot\alpha} = 4ma^2\dot\alpha,\qquad \frac{\partial\mathcal{L}}{\partial\alpha} = -2mga\alpha + ka^2(\beta-\alpha)\quad\Longrightarrow\quad 4ma^2\ddot\alpha + (2mga+ka^2)\alpha - ka^2\beta = 0$$

After dividing by $4ma^2$, all three equations (34), (36), (37) coincide exactly with (24)–(26). The two methods confirm each other.

Simulink models

The paper builds two models with three integrator chains each:

Nonlinear (equations 13, 17, 20): differences of the angles are formed with sum blocks, passed through sin blocks and gains $k/4m$. Gravity enters through a sin block and gain $g/2a$.
Linear (equations 24–26): only gains, $g/2a + k/4m$ for the outer pendulums and $g/2a + k/2m$ for the middle one, with coupling gains $k/4m$.

The mass matrix is diagonal, so there are no algebraic loops. The results are sent to the MATLAB workspace and plotted by a script that asks for the parameters and lets the user compare the two models.

Parameters used: $m = 1$ kg, $a = 1.5$ m, $g = 9.80665$ m/s², $k = 3.25$ N/m, with initial conditions $\alpha_0 = \beta_0 = 0$ and $\gamma_0 = \pi/10$ (18°), all velocities zero.

Simulation results from the paper

Over 20 s, the pendulum released first (γ, starting at 0.314 rad) gradually loses amplitude. The other two pendulums pick it up: α reaches about 0.29 rad at $t\approx10.7$ s, when γ and β are almost at rest. The energy then flows back. The paper plots the nonlinear and linear responses one above the other and concludes that there are no deviations between them, so the approximation is valid for small angles.

The simulator reproduces these plots: for example $\alpha(10.7\,\text{s}) = 0.286$ rad.

Review and normal-mode analysis

The derivations are correctBoth the torque derivation (13, 17, 20) and the Lagrangian derivation (34, 36, 37) are correct, and the linear equations agree. Only small slips appear: $d\omega/dT$ instead of $d\omega/dt$, and "generalisani" spellings.

7.1 Three different "nonlinear" spring models

The two derivations use different spring models, and neither is the exact one:

Torque method: force $=k(\overrightarrow{AB}-\overrightarrow{O_1O_2})$, proportional to the change of the whole connecting vector. Potential energy: $ka^2\big(1-\cos(\alpha-\beta)\big)$, torque $a^2k\sin(\alpha-\beta)$.
Lagrange method: only the horizontal separation stretches the spring. Potential energy: $\tfrac12ka^2(\sin\beta-\sin\alpha)^2$.
A real spring of free length $d$ (the pivot spacing): the force acts along the spring with magnitude $k(|AB|-d)$.

All three have the same linearization, so for small angles they are equivalent. At larger angles they differ, and the simulator offers all of them. The paper's "nonlinear" model is the first one.

7.2 "No deviations" – true only approximately

At 18° the two models are close but not identical. The phase difference accumulates over time, and by 20 s the maximum difference is about 0.028 rad (around 10 % of the amplitude). At 45° the difference grows to about 0.43 rad and the curves look completely different. The conclusion should be stated quantitatively: the approximation error grows with amplitude and with time.

7.3 Normal modes explain the energy exchange

In matrix form the linear system is $\ddot{\boldsymbol\theta} + \left(\frac{g}{2a}\mathbf I + \frac{k}{4m}\mathbf P\right)\boldsymbol\theta = \mathbf 0$, where $\mathbf P = \begin{bmatrix}1&-1&0\\-1&2&-1\\0&-1&1\end{bmatrix}$ has eigenvalues 0, 1 and 3:

Mode	Shape [α, β, γ]	ω²	ω (paper's values)	What happens
1	[1, 1, 1]	$g/2a$	1.808 rad/s	all swing together, springs never stretch
2	[1, 0, −1]	$g/2a + k/4m$	2.020 rad/s	outer pendulums in opposition, middle at rest
3	[1, −2, 1]	$g/2a + 3k/4m$	2.389 rad/s	middle opposes the outer two

Releasing only γ excites all three modes. Their frequencies are close (the springs are weak compared with gravity), so the modes drift in and out of phase. That is the slow energy transfer from γ to α and back seen in the plots: a three-pendulum version of beats. Start the system in one mode shape and it oscillates at a single frequency forever, with no transfer.

Interactive simulator

The simulator solves the equations numerically (4th-order Runge–Kutta, 2 ms step). Angles are entered in degrees, and plots show radians as in the paper. For the "real spring" model the pivot spacing $d$ sets the free length of the springs.

Cases

Model

Compare with

System parameters

m [kg]

a (half length) [m]

k [N/m]

g [m/s²]

d (pivot spacing) [m]

Initial angles [°] (velocities 0)

α₀

β₀

γ₀

Duration [s]

Angles α, β, γ

Difference between the two models

Energy of each pendulum (selected model)

Kinetic + gravitational energy of each pendulum: shows energy moving from one pendulum to another.

Animation (selected model)

Results table (every 0.1 s)

Review questions

Why is the moment of inertia $4ma^2$?

The rods are massless and the whole mass $m$ is concentrated at distance $2a$ from the pivot: $I = m(2a)^2$.

Why is the Lagrangian method "shorter and easier" here?

It needs only two scalar energies; no vector cross products and no bookkeeping of force directions. Internal forces (the springs acting on both neighbours) are handled automatically through the potential energy.

In which mode do the springs never stretch, and what is its frequency?

Mode [1, 1, 1]: all pendulums swing together, so the spring ends move identically. Its frequency is that of a single pendulum, $\omega = \sqrt{g/2a}$, independent of $k$.

Why does the energy move from pendulum γ to pendulum α?

The initial state is a superposition of three modes with slightly different frequencies. As the phase differences between them change, the modes add up at α and cancel at γ, and later the other way round. The weaker the springs, the closer the frequencies and the slower the exchange.

Why can the three spring models give different nonlinear results?

They agree in the first-order terms (horizontal displacement $a\theta$) but differ in the higher-order terms: vertical displacement of the attachment points, the spring's tilt, and its free length. Those terms matter only at larger angles.

How would you state the validity of the small-angle model quantitatively?

Compare both models over the simulation time and report the maximum difference relative to the amplitude. For example: 18° gives about 10 % after 20 s, while 45° gives more than 50 %. The error grows with both amplitude and time, because the frequencies differ slightly and the phase error accumulates.

Mode	Shape [α, β, γ]	ω²	ω (paper's values)	What happens
1	[1, 1, 1]	\(g/2a\)	1.808 rad/s	all swing together, springs never stretch
2	[1, 0, −1]	\(g/2a + k/4m\)	2.020 rad/s	outer pendulums in opposition, middle at rest
3	[1, −2, 1]	\(g/2a + 3k/4m\)	2.389 rad/s	middle opposes the outer two

Explanation of Method