A hanging chain

I remember reading Feynman II back in my youth and being blown away by the calculus of variations chapter, and immediately trying to play with it in every single different context I could. This is one of them, slightly more complex than those in the gospel according to RPF. We’re going to go on a route around the world and end up in a perhaps unexpected place.

Suppose you’ve got a chain hanging under gravity. What shape does it take? Let’s parametrize the shape as $(𝑥 (𝑠), 𝑦 (𝑠))$ where $𝑠$ is arc length. Then its potential energy is

\int_{0}^{𝐿} 𝜌 𝑔 𝑦 (𝑠) d 𝑠

where $𝐿$ is total length, $𝜌$ is mass-per-unit-length and $𝑔$ is, of course, gravity. When the chain is in steady state this will be at a minimum.

But, I hear you say, this clearly has no minimum because $𝑦$ can be made arbitrarily negative. There is a trick — remember I said that $𝑠$ is arc length. We need to include that constraint in our minimization; what we actually minimize is

ℒ︀ = \int_{0}^{𝐿} 𝜌 𝑔 𝑦 (𝑠) + \frac{𝜆 (𝑠)}{2} ({(\frac{d 𝑥}{d 𝑠})}^{2} + {(\frac{d 𝑦}{d 𝑠})}^{2} - 1) d 𝑠

where $𝜆 (𝑠)$ is a Lagrange multiplier to enforce the arc length constraint¹. It’ll actually turn out to be the tension in the chain.

Now the calculus of variations trick. What we’re going to do is to say that $ℒ︀$ is stationary with respect to first-order variations in $𝑥 (𝑠)$ , $𝑦 (𝑠)$ , and $𝜆 (𝑠)$ . This is exactly the same “stationary to first order” trick that works for minimization of functions, here played on this integral. Let’s set $𝑥 (𝑠) \mapsto 𝑥 (𝑠) + 𝛿 𝑥 (𝑠)$ , $𝑦 (𝑠) \mapsto 𝑦 (𝑠) + 𝛿 𝑦 (𝑠)$ , and $𝜆 (𝑠) \to 𝜆 (𝑠) + 𝛿 𝜆 (𝑠)$ and catch the first order terms. We get

𝛿 ℒ︀ = \int_{0}^{𝐿} 𝜌 𝑔 𝛿 𝑦 (𝑠) + \frac{𝛿 𝜆 (𝑠)}{2} ({(\frac{d 𝑥}{d 𝑠})}^{2} + {(\frac{d 𝑦}{d 𝑠})}^{2} - 1) + 𝜆 (𝑠) (\frac{d 𝑥}{d 𝑠} \frac{d 𝛿 𝑥}{d 𝑠} + \frac{d 𝑦}{d 𝑠} \frac{d 𝛿 𝑦}{d 𝑠}) d 𝑠,

and a bit of staring reveals that we’ve got $\frac{d 𝛿 𝑥}{d 𝑠}$ and $\frac{d 𝛿 𝑦}{d 𝑠}$ terms that we don’t know what to do with. So we need to get rid of them. Fortunately we remember that we can integrate by parts to get

𝛿 ℒ︀ = \int_{0}^{𝐿} (𝜌 𝑔 - \frac{d}{d 𝑠} (𝜆 (𝑠) \frac{d 𝑦}{d 𝑠})) 𝛿 𝑦 (𝑠) - (\frac{d}{d 𝑠} (𝜆 (𝑠) \frac{d 𝑥}{d 𝑠})) 𝛿 𝑥 (𝑠) + ({(\frac{d 𝑥}{d 𝑠})}^{2} + {(\frac{d 𝑦}{d 𝑠})}^{2} - 1) \frac{𝛿 𝜆 (𝑠)}{2} d 𝑠 .

The endpoint terms go away because we require $𝛿 𝑥 (0) = 𝛿 𝑦 (0) = 𝛿 𝑥 (𝐿) = 𝛿 𝑦 (𝐿) = 0$ — the chain is fixed at its start and end and our variations $𝛿 𝑥$ and $𝛿 𝑦$ must respect that.

Differential equations for fun and profit

Our variational terms are otherwise arbitrary, so if we want first-order variations to be zero (i.e. stationary energy) we get

\begin{aligned} \frac{d}{d 𝑠} (𝜆 (𝑠) \frac{d 𝑦}{d 𝑠}) & = 𝜌 𝑔, \\ \frac{d}{d 𝑠} (𝜆 (𝑠) \frac{d 𝑥}{d 𝑠}) & = 0, \\ {(\frac{d 𝑥}{d 𝑠})}^{2} + {(\frac{d 𝑦}{d 𝑠})}^{2} & = 1 . \end{aligned}

We’ll impose the conditions $𝑦 (0) = 𝑦 (𝐿) = 0$ , $𝑥 (0) = - 𝑎$ , $𝑥 (𝐿) = 𝑎$ , so that our string hangs between $(- 𝑎, 0)$ and $(𝑎, 0)$ . Integrating the first two equations we get

\begin{aligned} 𝜆 (𝑠) \frac{d 𝑥}{d 𝑠} & = 𝜌 𝑔 𝑝, \\ 𝜆 (𝑠) \frac{d 𝑦}{d 𝑠} & = 𝜌 𝑔 (𝑠 - 𝑠_{*}) . \end{aligned}

in which $𝑝$ and $𝑠_{*}$ are constants of integration². Now we substitute into the arc length constraint to get

{𝜆 (𝑠)}^{2} = 𝜌^{2} 𝑔^{2} (𝑝^{2} + {(𝑠 - 𝑠_{*})}^{2})

and therefore³

\frac{d 𝑦}{d 𝑠} = \frac{𝑠 - 𝑠_{*}}{\sqrt{𝑝^{2} + {(𝑠 - 𝑠_{*})}^{2}}},

giving

𝑦 (𝑠) = 𝑞 + \sqrt{𝑝^{2} + {(𝑠 - 𝑠_{*})}^{2}} .

We are almost home and dry:

\frac{d 𝑥}{d 𝑠} = \frac{1}{\sqrt{1 + 𝑝^{- 2} {(𝑠 - 𝑠_{*})}^{2}}} .

Make the change of variable $𝑠 = 𝑠_{*} + 𝑝 sinh 𝜃$ to get $𝑥 = 𝑝 (𝜃 - 𝜃_{0})$ , which we rewrite as

𝑦 = 𝑞 + 𝑝 cosh (\frac{𝑥}{𝑝} + 𝜃_{0}) .

Continuing on our merry way we find $𝜃_{0} = 0$ and

𝑦 = 𝑝 (cosh \frac{𝑥}{𝑝} - cosh \frac{𝑎}{𝑝}) .

We get $𝑝$ by solving $𝐿 = 2 𝑝 sinh \frac{𝑎}{𝑝}$ ; the simplest way to do this is to let $𝑝 = \frac{𝑎}{𝜉}$ and solve $\frac{𝐿}{2 𝑎} = \frac{sinh 𝜉}{𝜉}$ . Reassuringly, we see there is no solution if $𝐿 < 2 𝑎$ . A perhaps slightly surprising result is that the shape does not depend on $𝜌$ or $𝑔$ — all heavy chains hang the same way. (This is in fact not surprising: by redefining $𝜆 (𝑠)$ we can factor $𝜌 𝑔$ out of our Lagrangian.)

(Some years ago those nice guys Vella, Metcalfe, and Whittaker did something very much like this for floating rafts.)

Where this goes

Honestly, lots of places. One route ends up in Serious Numerical Analysis — the finite element methods that are the standard PDE solvers for solid mechanics and the Sobolev spaces that are the mathematical theory thereof⁴. Another route ends up in the very beautiful Lagrangian and Hamiltonian formulations of physics.

And one trivial and irrelevant byway that second route goes down is symplectic geometry, where you learn that the iteration

\begin{matrix} 𝑝_{𝑛 + \frac{1}{2}} = 𝑝_{𝑛} - {\frac{𝜕 𝑈}{𝜕 𝑞} |}_{𝑞_{𝑛}} \frac{Δ 𝑡}{2} \\ 𝑞_{𝑛 + 1} = 𝑞_{𝑛} + 𝑝_{𝑛 + \frac{1}{2}} Δ 𝑡 \\ 𝑝_{𝑛 + 1} = 𝑝_{𝑛 + \frac{1}{2}} - {\frac{𝜕 𝑈}{𝜕 𝑞} |}_{𝑞_{𝑛 + 1}} \frac{Δ 𝑡}{2} \end{matrix}

exactly preserves the 2-form $d 𝑝_{𝑛} \land d 𝑞_{𝑛}$ . And that is the reason that Hamiltonian Monte Carlo methods are so effective, linking us neatly back to my usual beat⁵. Maths… gotta know it all.

¹If I was a real man I would use $𝜆 (𝑠) (\sqrt{{(\frac{d 𝑥}{d 𝑠})}^{2} + {(\frac{d 𝑦}{d 𝑠})}^{2}} - 1)$ , but we will in fact get the same differential equations back from our version.
²I have been cunning and scaled the constants appropriately because I know how this is going to go.
³There is a choice of sign here: our optimization has two solutions — a maximum in which the chain goes up and a minimum in which the chain goes down. We want the minimum and have chosen the positive square root, and, slight subtlety, this means that $𝑝 > 0$ .
⁴If you try solving our problem this way, note that our Lagrangian integral needs $𝑥 (𝑠)$ and $𝑦 (𝑠)$ to have square integrable derivatives, and $𝜆 (𝑠)$ just has to be integrable. The simplest way of doing this is to let $𝑥 (𝑠)$ and $𝑦 (𝑠)$ be piecewise linear, and $𝜆 (𝑠)$ be piecewise constant. You’ll find (once you factor $𝜌 𝑔$ out of $𝜆$ )
$\begin{aligned} 0 & = \frac{𝜆_{𝑖 - \frac{1}{2}}}{𝛿_{𝑖 - \frac{1}{2}}} (𝑥_{𝑖} - 𝑥_{𝑖 - 1}) + \frac{𝜆_{𝑖 + \frac{1}{2}}}{𝛿_{𝑖 + \frac{1}{2}}} (𝑥_{𝑖} - 𝑥_{𝑖 + 1}), \\ 0 & = \frac{1}{2} 𝛿_{𝑖 - \frac{1}{2}} + \frac{1}{2} 𝛿_{𝑖 + \frac{1}{2}} + \frac{𝜆_{𝑖 - \frac{1}{2}}}{𝛿_{𝑖 - \frac{1}{2}}} (𝑦_{𝑖} - 𝑦_{𝑖 - 1}) + \frac{𝜆_{𝑖 + \frac{1}{2}}}{𝛿_{𝑖 + \frac{1}{2}}} (𝑦_{𝑖} - 𝑦_{𝑖 + 1}), \\ 1 & = \frac{{(𝑥_{𝑖 + 1} - 𝑥_{𝑖})}^{2}}{𝛿_{𝑖 + \frac{1}{2}}^{2}} + \frac{{(𝑦_{𝑖 + 1} - 𝑦_{𝑖})}^{2}}{𝛿_{𝑖 + \frac{1}{2}}^{2}} . \end{aligned}$
which you solve on $0 = 𝑠_{0} < 𝑠_{1} < \dots < 𝑠_{𝑁 - 1} < 𝑠_{𝑁} = 𝐿$ with $𝑥_{0} = - 𝑎$ , $𝑥_{𝑁} = 𝑎$ , $𝑦_{0} = 𝑦_{𝑁} = 0$ probably by using Newton iteration from some initial guess (the linear algebra is sparse).
⁵When the AI folk start muttering about “hidden structure on the latent manifold” I feel a strong urge to describe HMC as “iterated stochastic symplectic transformations on the cotangent bundle of the latent manifold”.