## Saturday, 29 June 2013

This entry is quite long; skip to the end to look at a piece of string if you don't feel like doing any maths. 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, it is completely equivalent to Newton's equations. 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.

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 = \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. Below is a simulation of N coupled pendulua for $N = 1 \ldots 12$. Presumably you could let $N \rightarrow \infty$ and get a continuous string, I'll leave that to you.