This entry is quite mathematical, remember the advice from last time about skipping equations, the text should be enough. The point of this is to introduce a useful object called the Lagrangian. The Lagrangian is the difference between kinetic and potential energy terms.
\[
L = T – V
\]
All of classical mechanics can be done using the Lagrangian and the equations of motion it obeys, which are completely equivalent to Newton’s equations. The animation below works by clicking somewhere in the box to start the pendulum. The slider for $N$ controls the number of points and the $\beta$ slider controls air-resistance.

If you have never seen a Lagrangian before you will have to take my word that it satisfies the equation
\[
\frac{d}{dt} \frac{\partial L}{\partial \dot{\theta}_i } = \frac{\partial L}{\partial \theta_i }
\]
where the dot denotes a time derivative and the $\theta_i$ are the co-ordinates relevant for the problem. We’re going to use pendulums which is why I chose the symbol $\theta$. Deriving this takes a while and involves concepts like D’Alembert’s principle principle and variational calculus.
Writing the equations of classical mechanics in this way is the beginning of a long story culminating in the principle of least action and the Hamiltonian. This in turn is how you start to write down quantum mechanics and quantum field theory. In fact the principle of least action which seems a bit mysterious in classical mechanics finds an explaination in quantum theory.

Now, following here, we start a quite involved derivation of the equations of motion for N coupled pendula. If you never studied calculus before I absolve you of guilt for not reading further. Even if you have it is still difficult, this is generally something 3^{rd} or 4^{th} year physics undergraduates might do. If you happen to be a 3^{rd} or 4^{th} year physics undergraduate with some free time, let’s begin!

Accepting Lagrange’s equations we can expand out the $\frac{d}{dt}$ to get,
\[
\sum_j \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \theta_j }\dot{\theta}_j
+ \sum_j \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \dot{\theta}_j }\ddot{\theta}_j
= \frac{\partial L}{\partial \theta_i }.
\]
The second partial derivatives have two indices, $i$ and $j$ and act on the one index object $\theta_j$ in the same way a matrix multiplies a vector. We define the two matrices
\[
M_{ij} = \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \dot{\theta}_j }
\]
the “mass” (for zero potential this would just be a diagonal matrix with the particle masses on the diagonal) and
\[
F_{ij} = \frac{\partial^2 L}{\partial \dot{\theta}_i \partial \theta_j }
\]
The equation is now,
\[
F \dot{\theta} + M \ddot{ \theta } = \frac{\partial L}{\partial \theta}
\]
I have dropped the indices and you should think of $\theta$ as a vector.

As in the previous post for numerical analysis we prefer first order equations, even if there are more of them. We define $\omega = \dot{\theta} $ and re-write the equation as,
\[
\frac{d}{dt} \left( \begin{array}{c}
\theta \\
\omega
\end{array} \right) =
\left( \begin{array}{cc}
0 & I \\
0 & -M^{-1}F
\end{array} \right)
\left( \begin{array}{c}
\theta \\
\omega
\end{array} \right) +
\left( \begin{array}{c}
0 \\
M^{-1} \frac{\partial L}{\partial \theta}
\end{array} \right)
\]
As bad as it looks this is now in a form amenable to numerical analysis. We just have to calculate, $M$, $F$ and $\frac{\partial L}{\partial \theta}$. Let’s get to it!

We’ll make things a bit simpler for ourselves by making all the pendulums of equal length and all the pendulum bobs of equal mass. The following diagram summarises the situation for two pendulums.

From the diagram we can see (after staring at it for a while) that
\[
\begin{array}{cc}
x_1 = l \sin \theta_1 & y_1 = L \cos \theta_1 \\
x_2 = l (\sin \theta_1+\sin \theta_2) & y_2 = l (\cos \theta_1 + \cos \theta_2) \\
\vdots & \vdots
\end{array}
\]
We will use the abbreviations
\[
\begin{array}{c}
c_i = \cos \theta_i \\
s_i = \sin \theta_i
\end{array}
\]
The potential of any of the masses is $-mgh$ where $h$ is the height. This gives for all $N$ pendula a potential,
\[
\begin{array}{l}
V = -mg\left[ y_1 + y_2 + \ldots + y_N \right] \\
V = -mgl\left[ c_1 + (c_1 + c_2) + \ldots + (c_1 + \ldots + c_N) \right] \\
V = -mgl\left[ N c_1 + (N-1) c_2 + \ldots + c_N \right]
\end{array}
\]
The kinetic energy is
\[
T = \frac{m}{2} \left[ \dot{x}_1^2 + \dot{y}_1^2 + \dot{x}_2^2 + \dot{y}_2^2 + \ldots + \dot{x}_N^2 + \dot{y}_N^2 \right].
\]
From the expressions for $x_i$ and $y_i$ we obtain,
\[
\begin{array}{cc}
\dot{x}_1 = l \cos \theta_1 \dot{\theta}_1 & \dot{y}_1 = l \sin \theta_1 \dot{\theta}_1\\
\dot{x}_2 = l (\cos \theta_1 \dot{\theta}_1 + \cos \theta_2 \dot{\theta}_2 ) & \dot{y}_2 = l (\sin \theta_1 \dot{\theta}_1 + \sin \theta_2 \dot{\theta}_2 ) \\
\vdots & \vdots
\end{array}
\]
and the kinetic energy,
\[
\begin{array}{l}
T = \frac{m}{2} \left[ \dot{x}_1^2 + \dot{y}_1^2 + \dot{x}_2^2 + \dot{y}_2^2 + \ldots + \dot{x}_N^2 + \dot{y}_N^2 \right].\\
T = \frac{ml^2}{2} \left[ (c_1 \omega_1)^2 + (c_1 \omega_1 + c_2 \omega_2 )^2 + \ldots + (c_1 \omega_1 + c_2 \omega_2 + \ldots c_N \omega_N ) + c \leftrightarrow s \right].
\end{array}
\]

This is a bit complicated so we will use matrix notation to keep track of things. With
\[
J_{ab} =
\begin{array}{ll}
(N+1 – b) & \quad a < b \\
(N+1 - a) & \quad a \geq b
\end{array}
\]
a symmetric matrix (write it out) we have,
\[
T = \frac{ml^2}{2} \sum_{a = 1}^{N} \sum_{b = 1}^{N} J_{ab} c_a c_b \omega_a \omega_b+ c \leftrightarrow s = \frac{ml^2}{2} \sum_{a = 1}^{N} \sum_{b = 1}^{N} J_{ab} c_{ab} \omega_a \omega_b
\]
where
\[
\begin{array}{c}
c_{ab} = \cos (\theta_a - \theta_b) = c_{ba} \\
s_{ab} = \sin (\theta_a - \theta_b) = -s_{ba}
\end{array}
\]

Let’s start by calculating $M$ and $F$.
\[
\frac{ \partial T}{\partial \omega_j } = ml^2 \sum_{a = 1}^{N} \omega_a J_{aj} c_{aj}
\]
Then
\[
M_{ij} = \frac{ \partial^2 T}{\partial \omega_i \partial \omega_j } = ml^2 J_{ij} c_{ij}
\]
and using the symetry of mixed partial derivatives
\[
F_{ij} = \frac{ \partial^2 T}{\partial \theta_j \partial \omega_i } = ml^2 \left[ J_{ij} \omega_j s_{ij} – \delta_{ij} \sum_{a = 1}^{N} \omega_a J_{ai} s_{ia} \right]
\]
Now we need $\frac{\partial L}{\partial \theta_i}$. We have
\[
\frac{\partial V}{\partial \theta_i} = mgl J_{1i} s_i
\]
and
\[
\frac{\partial T}{\partial \theta_i} = ml^2 \sum_{a=1}^N J_{ai} s_{ai} \omega_a \omega_i
\]
Which gives $\frac{\partial L}{\partial \theta_i}$. Phew!

There actually is one more term. Adding damping in the Lagrangian framework is a little bit tricky, a full explanation is here. This gives us an extra parameter $\beta$ which we use to control the amount of drag. You could let $N \rightarrow \infty$ and get a continuous string but I’ll leave that to you.