Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

 papasu on PHY450H1S. Relativistic Electr… papasu on Energy term of the Lorentz for… lidiodu on PHY450H1S. Relativistic Electr… lidiodu on PHY450H1S. Relativistic Electr… lidiodu on bivector form of Stokes t…

• 310,898

Archive for October, 2009

some rudimentary SQL

Posted by peeterjoot on October 28, 2009

Ironically, despite 10 years of working on a database, my work has been so far in the guts that I rarely have to do any SQL myself.  At one point when I’d been drafted for some some project management related work I did some fancy stuff with it making forms with python to compare regression results, but I’ve already forgotten all that.  Here’s a couple notes to myself of the very simplest stuff.

— create a table.  Note the backwards “prototype” compared to a C function;)

create table b (x int)

— how many rows?

select count(*) from b

— double table size

insert into b select * from b

–touch the whole thing to force IO and general activity in the database engine

update b set x = (x+1)

clearcase vs. /proc//fd/. clearcase within setview looses.

Posted by peeterjoot on October 27, 2009

Here’s a curious clash of virtual filesystems. Am trying to access my own processes’ /proc/-pid-/fd directory to investigate a file descriptor leak, and am unable to do so:

$ps -ef | grep db2sysc | grep peeter | grep -v grep | tail -1 peeterj 9318 9316 99 12:09 ? 01:03:14 db2sysc 0$  cd /proc/9318/fd
bash: cd: /proc/9318/fd: Permission denied
$cd /proc/9318$  ls
auxv  cpuset   environ  fd   mapped_base  mem   numa_maps  oom_score  seccomp  stat   status  wchan

I'd actually seen this before because we have code in our product that tries to access /proc/-pid-/stat stuff, and it doesn't work properly (sometimes and mysteriously).  Even odder, I can't even get at this as root
# ps -o pid -o ruid -o euid -o suid -o fsuid -o fname -a | grep 
21861     0     0     0     0 sh
# cd /proc/9318/fd
sh: cd: /proc/9318/fd: Permission denied
# cd /proc/9318
# ls
auxv  cpuset   environ  fd   mapped_base  mem   numa_maps  oom_score  seccomp  stat   status  wchan
# ls -l
total 0
dr-xr-xr-x   2 peeterj pdxdb2 0 2009-10-27 12:11 attr
-r--------   1 peeterj pdxdb2 0 2009-10-27 12:10 auxv


Something funny is happening in the kernel, since my session does appear to have sufficient root-ish behaviour (even the linux filesystem fsuid is set right). Turns out that this is some kind of clash between the clearcase version control virtual filesystem and the /proc virtual filesystem. When I am in my view, even as root:

# /usr/atria/bin/cleartool pwv
Working directory view: ** NONE **
Set view: peeterj_o26
#

I have no access to much of /proc/, but running as any old user when there is no trouble

$/usr/atria/bin/cleartool pwv Working directory view: ** NONE ** Set view: ** NONE **$  pwd
/proc/9318/fd

What a bizarre quirk! Glad to have this figured out … now back to the file descriptor leak.

A stupid coding error for today.

Posted by peeterjoot on October 27, 2009

I meant for the second if to be an else-if, and was scratching my head for a while trying to figure out how I got to the “impossible” UNEXPECTED codepath. Ugg.

inline sqlzRc MapReturnCode( int     nativeRc )
{
sqlzRc rc = SQLO_OK ;

if ( SUCCESS == nativeRc )
{
rc = SQLO_OK ;
}
if ( IsErrorFoo( nativeRc ) )
{
rc = SQLE_FOO ;
}
else if ( IsErrorBlah( nativeRc ) )
{
rc = SQLE_BLAH ;
}
...
else
{
rc = SQLE_UNEXPECTED_ERROR ;


spherical polar line element (continued).

Posted by peeterjoot on October 25, 2009

[Click here for a PDF of this sequence of posts with nicer formatting]

Line element.

The line element for the particle moving on a spherical surface can be calculated by calculating the derivative of the spherical polar unit vector $\hat{\mathbf{r}}$

\begin{aligned}\frac{d\hat{\mathbf{r}}}{dt} = \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \frac{d\phi}{dt} +\frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \frac{d\theta}{dt}\end{aligned} \quad\quad\quad(17)

than taking the magnitude of this vector. We can start either in coordinate form

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta \cos\phi + \mathbf{e}_2 \sin\theta \sin\phi\end{aligned} \quad\quad\quad(18)

or, instead do it the fun way, first grouping this into a complex exponential form. This factorization was done above, but starting over allows this to be done a bit more effectively for this particular problem. As above, writing $i = \mathbf{e}_1 \mathbf{e}_2$, the first factorization is

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta e^{i\phi} \end{aligned} \quad\quad\quad(19)

The unit vector $\boldsymbol{\rho} = \mathbf{e}_1 e^{i\phi}$ lies in the $x,y$ plane perpendicular to $\mathbf{e}_3$, so we can form the unit bivector $\mathbf{e}_3\boldsymbol{\rho}$ and further factor the unit vector terms into a $\cos + i \sin$ form

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta e^{i\phi} \\ &= \mathbf{e}_3 (\cos\theta + \mathbf{e}_3 \boldsymbol{\rho} \sin\theta) \\ \end{aligned}

This allows the spherical polar unit vector to be expressed in complex exponential form (really a vector-quaternion product)

\begin{aligned}\hat{\mathbf{r}} = \mathbf{e}_3 e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} = e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \mathbf{e}_3\end{aligned} \quad\quad\quad(20)

Now, calculating the unit vector velocity, we get

\begin{aligned}\frac{d\hat{\mathbf{r}}}{dt} &= \mathbf{e}_3 \mathbf{e}_3 \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \dot{\theta} + \mathbf{e}_1 \mathbf{e}_1 \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \\ &= \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \left(\dot{\theta} + e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho} \sin\theta e^{-i\phi} \mathbf{e}_2 \dot{\phi}\right) \\ &= \left( \dot{\theta} + \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \right) e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho}\end{aligned}

The last two lines above factor the $\boldsymbol{\rho}$ vector and the $e^{\mathbf{e}_3 \boldsymbol{\rho} \theta}$ quaternion to the left and to the right in preparation for squaring this to calculate the magnitude.

\begin{aligned}\left( \frac{d\hat{\mathbf{r}}}{dt} \right)^2&=\left\langle{{ \left( \frac{d\hat{\mathbf{r}}}{dt} \right)^2 }}\right\rangle \\ &=\left\langle{{ \left( \dot{\theta} + \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \right) \left(\dot{\theta} + e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho} \sin\theta e^{-i\phi} \mathbf{e}_2 \dot{\phi}\right) }}\right\rangle \\ &=\dot{\theta}^2 + \sin^2\theta \dot{\phi}^2+ \sin\theta \dot{\phi} \dot{\theta}\left\langle{{ \mathbf{e}_2 e^{i\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \rho \theta}+e^{-\mathbf{e}_3 \rho \theta} \boldsymbol{\rho} e^{-i\phi} \mathbf{e}_2}}\right\rangle \\ \end{aligned}

This last term ($\in \text{span} \{\rho\mathbf{e}_1, \rho\mathbf{e}_2, \mathbf{e}_1\mathbf{e}_3, \mathbf{e}_2\mathbf{e}_3\}$) has only grade two components, so the scalar part is zero. We are left with the line element

\begin{aligned}\left(\frac{d (r\hat{\mathbf{r}})}{dt}\right)^2 = r^2 \left( \dot{\theta}^2 + \sin^2\theta \dot{\phi}^2 \right)\end{aligned} \quad\quad\quad(21)

In retrospect, at least once one sees the answer, it seems obvious. Keeping $\theta$ constant the length increment moving in the plane is $ds = \sin\theta d\phi$, and keeping $\phi$ constant, we have $ds = d\theta$. Since these are perpendicular directions we can add the lengths using the Pythagorean theorem.

Line element using an angle and unit bivector parameterization.

Parameterizing using scalar angles is not the only approach that we can take to calculate the line element on the unit sphere. Proceding directly with a alternate polar representation, utilizing a unit bivector $j$, and scalar angle $\theta$ is

\begin{aligned}\mathbf{x} = r \mathbf{e}_3 e^{j\theta}\end{aligned} \quad\quad\quad(22)

For this product to be a vector $j$ must contain $\mathbf{e}_3$ as a factor ($j = \mathbf{e}_3 \wedge m$ for some vector m.) Setting $r = 1$ for now, the deriviate of $\mathbf{x}$ is

\begin{aligned}\dot{\mathbf{x}} &= \mathbf{e}_3 \frac{d}{dt} \left( \cos\theta + j \sin\theta \right) \\ &= \mathbf{e}_3 \dot{\theta} \left( -\sin\theta + j \cos\theta \right) + \mathbf{e}_3 \frac{d j}{dt} \sin\theta \\ &= \mathbf{e}_3 \dot{\theta} j \left( j \sin\theta + \cos\theta \right) + \mathbf{e}_3 \frac{d j}{dt} \sin\theta \\ \end{aligned}

This is

\begin{aligned}\dot{\mathbf{x}} = \mathbf{e}_3 \left( \frac{d\theta}{dt} j e^{j\theta} + \frac{d j}{dt} \sin\theta \right)\end{aligned} \quad\quad\quad(23)

Alternately, we can take derivatives of $\mathbf{x} = r e^{-j\theta} \mathbf{e}_3$, for

Or

\begin{aligned}\dot{\mathbf{x}} = -\left( \frac{d\theta}{dt} j e^{-j\theta} + \frac{d j}{dt} \sin\theta \right) \mathbf{e}_3\end{aligned} \quad\quad\quad(24)

Together with 23, the line element for position change on the unit sphere is then

\begin{aligned}\dot{\mathbf{x}}^2 &= \left\langle{{-\left( \frac{d\theta}{dt} j e^{-j\theta} + \frac{d j}{dt} \sin\theta \right) \mathbf{e}_3 \mathbf{e}_3 \left( \frac{d\theta}{dt} j e^{j\theta} + \frac{d j}{dt} \sin\theta \right) }}\right\rangle \\ &= \left\langle{{-\left( \frac{d\theta}{dt} j e^{-j\theta} + \frac{d j}{dt} \sin\theta \right) \left( \frac{d\theta}{dt} j e^{j\theta} + \frac{d j}{dt} \sin\theta \right) }}\right\rangle \\ &=\left(\frac{d\theta}{dt}\right)^2 - \left(\frac{d j}{dt}\right)^2 \sin^2\theta - \frac{d\theta}{dt} \sin\theta \left\langle{{ \frac{dj}{dt} j e^{j\theta} + j e^{-j\theta} \frac{dj}{dt} }}\right\rangle \\ \end{aligned}

Starting with cyclic reordering of the last term, we get zero

\begin{aligned}\left\langle{{ \frac{dj}{dt} j e^{j\theta} + j e^{-j\theta} \frac{dj}{dt} }}\right\rangle &=\left\langle{{ \frac{dj}{dt} j \left( e^{j\theta} + e^{-j\theta} \right) }}\right\rangle \\ &=\left\langle{{ \frac{dj}{dt} j 2 j \sin\theta }}\right\rangle \\ &=- 2 \sin\theta \frac{d}{dt} \underbrace{\left\langle{{ j }}\right\rangle}_{=0} \\ \end{aligned}

The line element (for constant $r$) is therefore

\begin{aligned}\dot{\mathbf{x}}^2 =r^2 \left( \dot{\theta}^2 - \left(\frac{dj}{dt}\right)^2 \sin^2 \theta \right)\end{aligned} \quad\quad\quad(25)

This is essentially the same result as we got starting with an explicit $r, \theta, \phi$. Repeating for comparision that was

\begin{aligned}\dot{\mathbf{x}}^2 = r^2 \left( \dot{\theta}^2 + \sin^2\theta \dot{\phi}^2 \right)\end{aligned} \quad\quad\quad(26)

The bivector that we have used this time encodes the orientation of the plane of rotation from the polar axis down to the position on the sphere corresponds to the angle $\phi$ in the
scalar parameterization. The negation in sign is expected due to the negative bivector square.

Also comparing to previous results it is notable that we can explicitly express this bivector in terms of the scalar angle if desired as

\begin{aligned}\boldsymbol{\rho} &= \mathbf{e}_1 e^{\mathbf{e}_1 \mathbf{e}_2 \phi} = \mathbf{e}_1 \cos\phi + \mathbf{e}_2 \sin\phi \\ j &= \mathbf{e}_3 \wedge \boldsymbol{\rho} = \mathbf{e}_3 \boldsymbol{\rho}\end{aligned} \quad\quad\quad(27)

The inverse mapping, expressing the scalar angle using the bivector representation is also possible, but not unique. The principle angle for that inverse mapping is

\begin{aligned}\phi &= -\mathbf{e}_1 \mathbf{e}_2 \ln(\mathbf{e}_1 \mathbf{e}_3 j) \end{aligned} \quad\quad\quad(29)

Allowing the magnitude to vary.

Writing a vector in polar form

\begin{aligned}\mathbf{x} = r \hat{\mathbf{r}}\end{aligned} \quad\quad\quad(30)

and also allowing $r$ to vary, we have

\begin{aligned}\left(\frac{d\mathbf{x}}{dt}\right)^2 &= \left( \frac{dr}{dt} \hat{\mathbf{r}} + r \frac{d\hat{\mathbf{r}}}{dt} \right)^2 \\ &= \left( \frac{dr}{dt} \right)^2 + r^2 \left( \frac{d\hat{\mathbf{r}}}{dt} \right)^2 + 2 r \frac{dr}{dt} \hat{\mathbf{r}} \dot{c} \frac{d\hat{\mathbf{r}}}{dt} \end{aligned}

The squared unit vector derivative was previously calculated to be

\begin{aligned}\left(\frac{d \hat{\mathbf{r}}}{dt}\right)^2 = \dot{\theta}^2 + \sin^2\theta \dot{\phi}^2 \end{aligned} \quad\quad\quad(31)

Picturing the geometry is enough to know that $\dot{\hat{\mathbf{r}}} \dot{c} \hat{\mathbf{r}} = 0$ since $\dot{\hat{\mathbf{r}}}$ is always tangential to the sphere. It should also be possible to algebraically show this, but without going through the effort we at least now know the general line element

\begin{aligned}\dot{\mathbf{x}}^2 = {\dot{r}}^2 + r^2 \left( \dot{\theta}^2 + \sin^2\theta \dot{\phi}^2 \right)\end{aligned} \quad\quad\quad(32)

Hamiltonian equations for double and multiple planar pendulums, and general quadratic form velocity dependence.on

Posted by peeterjoot on October 19, 2009

(NOTE: the process of converting from standalone latex to wordpress got messed up in a few places. This have been fixed manually, but if anything seems wrong, check first in the standalone PDF above.)

In the following section I started off with the goal of treating two connected pendulums moving in a plane. Even setting up the Hamiltonians for this turned out to be a bit messy, requiring a matrix inversion. Tackling the problem in the guise of using a more general quadratic form (which works for the two particle as well as $N$ particle cases) seemed like it would actually be simpler than using the specifics from the angular velocity dependence of the specific pendulum problem. Once the Hamiltonian equations were found in this form, an attempt to do the first order Taylor expansion as done for the single planar pendulum and the spherical pendulum was performed. This turned out to be a nasty mess and is seen to not be particularily illuminating. I didn’t know that is how it would turn out ahead of time since I had my fingers crossed for some sort of magic simplification once the final substuitions were made. If such a simplification is possible, the procedure to do so is not obvious.

Although the Hamiltonian equations for a spherical pendulum have been considered previously, for the double pendulum case it seems prudent to avoid temptation, and to first see what happens with a simpler first step, a planar double pendulum.

Setting up coordinates $x$ axis down, and $y$ axis to the left with $i = \hat{\mathbf{x}} \hat{\mathbf{y}}$ we have for the position of the first mass $m_1$, at angle ${\theta_1}$ and length $l_1$

\begin{aligned}z_1 = \hat{\mathbf{x}} l_1 e^{i{\theta_1}}\end{aligned} \quad\quad\quad(180)

If the second mass, dangling from this is at an angle ${\theta_2}$ from the $x$ axis, its position is

\begin{aligned}z_2 = z_1 + \hat{\mathbf{x}} l_2 e^{i{\theta_2}}\end{aligned} \quad\quad\quad(181)

We need the velocities, and their magnitudes. For $z_1$ this is

\begin{aligned}{\left\lvert{\dot{z}_1}\right\rvert}^2 = {l_1}^2 {\dot{\theta}_1}^2\end{aligned} \quad\quad\quad(182)

For the second mass

\begin{aligned}\dot{z}_2 = \hat{\mathbf{x}} i \left( l_1 {\dot{\theta}_1} e^{i{\theta_1}} + l_2 {\dot{\theta}_2} e^{i{\theta_2}} \right)\end{aligned} \quad\quad\quad(183)

Taking conjugates and multiplying out we have

\begin{aligned}{\left\lvert{\dot{z}_2}\right\rvert}^2 ={l_1}^2 {\dot{\theta}_1}^2+ 2 l_1 l_2 {\dot{\theta}_1} {\dot{\theta}_2} \cos({\theta_1} - {\theta_2})+{l_2}^2 {\dot{\theta}_2}^2\end{aligned} \quad\quad\quad(184)

That’s all that we need for the Kinetic terms in the Lagrangian. Now we need the height for the $m g h$ terms. If we set the reference point at the lowest point for the double pendulum system, the height of the first particle is

\begin{aligned}h_1 = l_2 + l_1 (1 -\cos{\theta_1})\end{aligned} \quad\quad\quad(185)

For the second particle, the distance from the horizontal is

\begin{aligned}d = l_1 \cos{\theta_1} + l_2 \cos{\theta_2}\end{aligned} \quad\quad\quad(186)

So the total distance from the reference point is

\begin{aligned}h_2 = l_1 (1 - \cos{\theta_1}) + l_2 (1 -\cos{\theta_2})\end{aligned} \quad\quad\quad(187)

We now have the Lagrangian

\begin{aligned}\mathcal{L}' & =\frac{1}{{2}} m_1 {l_1}^2 {\dot{\theta}_1}^2+ \frac{1}{{2}} m_2 \left({l_1}^2 {\dot{\theta}_1}^2+ 2 l_1 l_2 {\dot{\theta}_1} {\dot{\theta}_2} \cos({\theta_1} - {\theta_2})+{l_2}^2 {\dot{\theta}_2}^2\right) \\ & -m_1 g (l_2 + l_1 (1 -\cos{\theta_1}))-m_2 g ( l_1 (1 - \cos{\theta_1}) + l_2 (1 -\cos{\theta_2}) )\end{aligned} \quad\quad\quad(188)

Dropping constant terms (effectively choosing a difference reference point for the potential) and rearranging a bit, also writing $M = m_1 + m_2$, we have the simpler Lagrangian

\begin{aligned}\mathcal{L} & =\frac{1}{{2}} M {l_1}^2 {\dot{\theta}_1}^2 + \frac{1}{{2}} m_2 {l_2}^2 {\dot{\theta}_2}^2+ m_2 l_1 l_2 {\dot{\theta}_1} {\dot{\theta}_2} \cos({\theta_1} - {\theta_2})+M l_1 g \cos{\theta_1}+m_2 l_2 g \cos{\theta_2}\end{aligned} \quad\quad\quad(190)

The conjugate momenta that we need for the Hamiltonian are

\begin{aligned}P_{\theta_1} & =M {l_1}^2 {\dot{\theta}_1}+ m_2 l_1 l_2 {\dot{\theta}_2} \cos({\theta_1} - {\theta_2}) \\ P_{\theta_2} & =m_2 {l_2}^2 {\dot{\theta}_2}+ m_2 l_1 l_2 {\dot{\theta}_1} \cos({\theta_1} - {\theta_2})\end{aligned} \quad\quad\quad(191)

Unlike any of the other simpler Hamiltonian systems considered so far, the coupling between the velocities means that we have a system of equations that we must first invert before we can even express the Hamiltonian in terms of the respective momenta.

That is

\begin{aligned}\begin{bmatrix}P_{\theta_1} \\ P_{\theta_2}\end{bmatrix}=\begin{bmatrix}M {l_1}^2 & m_2 l_1 l_2 \cos({\theta_1} - {\theta_2}) \\ m_2 l_1 l_2 \cos({\theta_1} - {\theta_2}) & m_2 {l_2}^2\end{bmatrix}\begin{bmatrix}{\dot{\theta}_1} \\ {\dot{\theta}_2}\end{bmatrix}\end{aligned} \quad\quad\quad(193)

While this is easily invertible, doing so and attempting to substitute it back, results in an unholy mess (albeit perhaps one that can be simplified). Is there a better way? A possibly promising way is motivated by observing that this matrix, a function of the angular difference $\delta = {\theta_1} - {\theta_2}$, looks like it is something like a moment of inertia tensor. If we call this $\mathcal{I}$, and write

\begin{aligned}\boldsymbol{\Theta} & \equiv \begin{bmatrix}{\theta_1} \\ {\theta_2}\end{bmatrix}\end{aligned} \quad\quad\quad(194)

Then the relation between the conjugate momenta in vector form

\begin{aligned}\mathbf{p} & \equiv \begin{bmatrix}P_{\theta_1} \\ P_{\theta_2}\end{bmatrix}\end{aligned} \quad\quad\quad(195)

and the angular velocity vector can be written

\begin{aligned}\mathbf{p} & = \mathcal{I}(\delta)\dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(196)

Can we write the Lagrangian in terms of $\dot{\boldsymbol{\Theta}}$? The first Kinetic term is easy, just

\begin{aligned}\frac{1}{{2}} m_1 l^2 {\dot{\theta}_1}^2 = \frac{1}{{2}} m_1\dot{\boldsymbol{\Theta}}^\text{T}\begin{bmatrix}l_1^2 & 0 \\ 0 & 0\end{bmatrix}\dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(197)

For the second mass, going back to 183, we can write

\begin{aligned}\dot{z}_2 = \hat{\mathbf{x}} i\begin{bmatrix}l_1 e^{i{\theta_1}} & l_2 e^{i{\theta_2}}\end{bmatrix}\dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(198)

Writing $\mathbf{r}$ for this $1x2$ matrix, we can utilize the associative property for compatible sized matrices to rewrite the speed for the second particle in terms of a quadratic form

\begin{aligned}{\left\lvert{\dot{z}_2}\right\rvert}^2=(\mathbf{r} \dot{\boldsymbol{\Theta}}) (\bar{\mathbf{r}} \dot{\boldsymbol{\Theta}})=\dot{\boldsymbol{\Theta}}^\text{T} (\mathbf{r}^\text{T} \bar{\mathbf{r}}) \dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(199)

The Lagrangian kinetic can all now be grouped into a single quadratic form

\begin{aligned}Q \equiv m_1\begin{bmatrix}l_1 \\ 0\end{bmatrix}\begin{bmatrix}l_1 & 0\end{bmatrix}+m_2\begin{bmatrix}l_1 e^{i{\theta_1}} \\ l_2 e^{i{\theta_2}}\end{bmatrix}\begin{bmatrix}l_1 e^{-i{\theta_1}} & l_2 e^{-i{\theta_2}}\end{bmatrix}\end{aligned} \quad\quad\quad(200)

\begin{aligned}\mathcal{L} =\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\text{T} Q \dot{\boldsymbol{\Theta}}+M l_1 g \cos{\theta_1}+m_2 l_2 g \cos{\theta_2}\end{aligned} \quad\quad\quad(201)

It is also clear that this generalize easily to multiple connected pendulums, as follows

\begin{aligned}K & = \frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\text{T} \sum_k m_k Q_k \dot{\boldsymbol{\Theta}} \\ Q_k & ={\begin{bmatrix}l_r l_s e^{i(\theta_r - \theta_s)}\end{bmatrix}}_{r,s \le k} \\ \phi & = - g \sum_k l_k \cos\theta_k \sum_{j=k}^N m_j \\ \mathcal{L} & = K - \phi\end{aligned} \quad\quad\quad(202)

In the expression for $Q_k$ above, it is implied that the matrix is zero for any indexes $r,s > k$, so it would perhaps be better to write explicitly

\begin{aligned}Q = \sum_k m_k Q_k={\begin{bmatrix}\sum_{j=\max(r,s)}^N m_j l_r l_s e^{i(\theta_r - \theta_s)}\end{bmatrix}}_{r,s}\end{aligned} \quad\quad\quad(206)

Returning to the problem, it is convenient and sufficient in many cases to only discuss the representative double pendulum case. For that we can calculate the conjugate momenta from 201 directly

\begin{aligned}P_{\theta_1} & =\frac{\partial {}}{\partial {{\dot{\theta}_1}}}\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\text{T} Q \dot{\boldsymbol{\Theta}} \\ & =\frac{\partial {}}{\partial {{\dot{\theta}_1}}}\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\text{T} Q\begin{bmatrix}1 \\ 0\end{bmatrix}+\frac{1}{{2}}\begin{bmatrix}1 & 0\end{bmatrix}Q \dot{\boldsymbol{\Theta}} \\ & =\begin{bmatrix}1 & 0\end{bmatrix}\left(\frac{1}{{2}}(Q + Q^\text{T})\right) \dot{\boldsymbol{\Theta}} \\ \end{aligned}

Similarly the ${\theta_2}$ conjugate momentum is

\begin{aligned}P_{\theta_2}=\begin{bmatrix}0 & 1\end{bmatrix}\left(\frac{1}{{2}}(Q + Q^\text{T})\right) \dot{\boldsymbol{\Theta}} \\ \end{aligned}

Putting both together, it is straightforward to verify that this recovers 193, which can now be written

\begin{aligned}\mathbf{p} & =\frac{1}{{2}}(Q + Q^\text{T}) \dot{\boldsymbol{\Theta}} = \mathcal{I} \dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(207)

Observing that $\mathcal{I} = \mathcal{I}^\text{T}$, and thus $(\mathcal{I}^\text{T})^{-1} = \mathcal{I}^{-1}$, we now have everything required to express the Hamiltonian in terms of the conjugate momenta

\begin{aligned}H = \mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1} \right) \mathbf{p} - M g l_1 \cos{\theta_1} - m_2 l_2 g \cos{\theta_2}\end{aligned} \quad\quad\quad(208)

This is now in a convenient form to calculate the first set of Hamiltonian equations.

\begin{aligned}\dot{\theta}_k & =\frac{\partial {H}}{\partial {P_{\theta_k}}} \\ & =\frac{\partial {\mathbf{p}^\text{T}}}{\partial {P_{\theta_k}}} \frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1} \mathbf{p}+ \mathbf{p}^\text{T} \frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1} \frac{\partial {\mathbf{p}^\text{T}}}{\partial {P_{\theta_k}}} \\ & ={\begin{bmatrix}\delta_{kj}\end{bmatrix}}_j\frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1} \mathbf{p}+ \mathbf{p}^\text{T} \frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1}{\begin{bmatrix}\delta_{ik}\end{bmatrix}}_i \\ & ={\begin{bmatrix}\delta_{kj}\end{bmatrix}}_j\mathcal{I}^{-1} \underbrace{\frac{1}{{2}}(Q + Q^\text{T})}_{\mathcal{I}} \mathcal{I}^{-1} \mathbf{p} \\ & ={\begin{bmatrix}\delta_{kj}\end{bmatrix}}_j\mathcal{I}^{-1} \mathbf{p} \\ \end{aligned}

So, when the velocity dependence is a quadratic form as identified in 200, the first half of the Hamiltonian equations in vector form are just

\begin{aligned}\dot{\boldsymbol{\Theta}} & ={\begin{bmatrix}\frac{\partial {}}{\partial {P_{\theta_1}}} \cdots \frac{\partial {}}{\partial {P_{\theta_N}}}\end{bmatrix}}^\text{T} H=\mathcal{I}^{-1} \mathbf{p}\end{aligned} \quad\quad\quad(209)

This is exactly the relation we used in the first place to re-express the Lagrangian in terms of the conjugate momenta in preparation for this calculation. The remaining Hamiltonian equations are trickier, and what we now want to calculate. Without specific reference to the pendulum problem, lets do this calculation for the general Hamiltonian for a non-velocity dependent potential. That is

\begin{aligned}H = \mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} Q \mathcal{I}^{-1} \right) \mathbf{p} + \phi(\boldsymbol{\Theta})\end{aligned} \quad\quad\quad(210)

The remaining Hamiltonian equations are ${\partial {H}}/{\partial {\theta_a}} = -\dot{P}_{\theta_a}$, and the tricky part of evaluating this is going to all reside in the Kinetic term. Diving right in this is

\begin{aligned}\frac{\partial {K}}{\partial {\theta_a}} & = \mathbf{p}^\text{T} \left( \frac{1}{{2}} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} Q \mathcal{I}^{-1} \right) \mathbf{p} +\mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} \frac{\partial {Q}}{\partial {\theta_a}} \mathcal{I}^{-1} \right) \mathbf{p} +\mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} Q \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} \right) \mathbf{p} \\ & = \mathbf{p}^\text{T} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} \underbrace{\frac{1}{{2}}(Q + Q^\text{T})}_{=\mathcal{I}} \mathcal{I}^{-1} \mathbf{p} +\mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} \frac{\partial {Q}}{\partial {\theta_a}} \mathcal{I}^{-1} \right) \mathbf{p} \\ & = \mathbf{p}^\text{T} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} \mathbf{p} +\mathbf{p}^\text{T} \left( \frac{1}{{2}} \mathcal{I}^{-1} \frac{\partial {Q}}{\partial {\theta_a}} \mathcal{I}^{-1} \right) \mathbf{p} \\ \end{aligned}

For the two particle case we can expand this inverse easily enough, and then take derivatives to evaluate this, but this is messier and intractable for the general case. We can however, calculate the derivative of the identity matrix using the standard trick from rigid body mechanics

\begin{aligned}0 & = \frac{\partial {I}}{\partial {\theta_a}} \\ & = \frac{\partial {(\mathcal{I} \mathcal{I}^{-1})}}{\partial {\theta_a}} \\ & = \frac{\partial {\mathcal{I}}}{\partial {\theta_a}} \mathcal{I}^{-1}+\mathcal{I} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} \\ \end{aligned}

Thus the derivative of the inverse (moment of inertia?) matrix is

\begin{aligned}\frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_a}} & = -\mathcal{I}^{-1}\frac{\partial {\mathcal{I}}}{\partial {\theta_a}} \mathcal{I}^{-1} \\ & = -\mathcal{I}^{-1}\frac{1}{{2}}\left( \frac{\partial {Q}}{\partial {\theta_a}} + \frac{\partial {Q^\text{T}}}{\partial {\theta_a}} \right) \mathcal{I}^{-1}\end{aligned}

This gives us for the Hamiltonian equation

\begin{aligned}\frac{\partial {H}}{\partial {\theta_a}} & = -\frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_a}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} + \frac{\partial {\phi}}{\partial {\theta_a}} \end{aligned} \quad\quad\quad(211)

If we introduce a phase space position gradients

\begin{aligned}\nabla & \equiv {\begin{bmatrix}\frac{\partial {}}{\partial {\theta_1}} \cdots \frac{\partial {}}{\partial {\theta_N}}\end{bmatrix}}^\text{T} \\ \end{aligned} \quad\quad\quad(212)

Then for the second half of the Hamiltonian equations we have the vector form

\begin{aligned}-\nabla H = \dot{\mathbf{p}} = {\begin{bmatrix}\frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \end{bmatrix}}_r- \nabla \phi\end{aligned} \quad\quad\quad(214)

The complete set of Hamiltonian equations for 210, in block matrix form, describing all the phase space change of the system is therefore

\begin{aligned}\frac{d}{dt} \begin{bmatrix}\mathbf{p} \\ \boldsymbol{\Theta}\end{bmatrix}=\begin{bmatrix}{\begin{bmatrix}\frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \end{bmatrix}}_r- \nabla \phi \\ \mathcal{I}^{-1} \mathbf{p} \\ \end{bmatrix}=\begin{bmatrix}{\begin{bmatrix}\frac{1}{{2}} \dot{\boldsymbol{\Theta}} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_r- \nabla \phi \\ \dot{\boldsymbol{\Theta}}\end{bmatrix}\end{aligned} \quad\quad\quad(215)

This is a very general relation, much more so than required for the original two particle problem. We have the same non-linearity that prevents this from being easily solved. If we want a linear expansion around a phase space point to find an approximate first order solution, we can get that applying the chain rule, calculating all the ${\partial {}}/{\partial {\theta_k}}$, and ${\partial {}}/{\partial {P_{\theta_k}}}$ derivatives of the top $N$ rows of this matrix.

If we write

\begin{aligned}\mathbf{z} \equiv \begin{bmatrix}\mathbf{p} \\ \boldsymbol{\Theta}\end{bmatrix}-{{\begin{bmatrix}\mathbf{p} \\ \boldsymbol{\Theta}\end{bmatrix}}}_{{t=0}}\end{aligned} \quad\quad\quad(216)

and the Hamiltonian equations 215 as

\begin{aligned}\frac{d}{dt}\begin{bmatrix}\mathbf{p} \\ \boldsymbol{\Theta}\end{bmatrix}= A(\mathbf{p}, \boldsymbol{\Theta})\end{aligned} \quad\quad\quad(217)

Then the linearization, without simplifying or making explicit yet is

\begin{aligned}\dot{\mathbf{z}} \approx{{\begin{bmatrix}{\begin{bmatrix}\frac{1}{{2}} \dot{\boldsymbol{\Theta}} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_r- \nabla \phi \\ \dot{\boldsymbol{\Theta}}\end{bmatrix}}}_{{t=0}}+{\left.{{\begin{bmatrix}\frac{\partial {A}}{\partial {P_{\theta_1}}} \cdots \frac{\partial {A}}{\partial {P_{\theta_N}}} \cdots \frac{\partial {A}}{\partial {\theta_1}} \cdots \frac{\partial {A}}{\partial {\theta_N}} \end{bmatrix}}}\right\vert}_{{t=0}} \mathbf{z}\end{aligned} \quad\quad\quad(218)

For brevity the constant term evaluated at $t=0$ is expressed in terms of the original angular velocity vector from our Lagrangian. The task is now evaluating the derivatives in the first order term of this Taylor series. Let’s do these one at a time and then reassemble all the results afterward.

So that we can discuss just the first order terms lets write $\Delta$ for the matrix of first order derivatives in our Taylor expansion, as in

\begin{aligned}f(\mathbf{p}, \boldsymbol{\Theta}) = {\left.{{f(\mathbf{p}, \boldsymbol{\Theta})}}\right\vert}_{{0}} + {\left.{{\Delta f}}\right\vert}_{{0}} \mathbf{z} + \cdots \end{aligned} \quad\quad\quad(219)

First, lets do the potential gradient.

\begin{aligned} \Delta (\nabla \phi) = \begin{bmatrix} 0 & {\begin{bmatrix} \frac{\partial^2 \phi}{ \partial \theta_r \partial \theta_c } \end{bmatrix}}_{r,c} \end{bmatrix} \end{aligned} \quad\quad\quad(220)

Next in terms of complexity is the first order term of $\dot{\boldsymbol{\Theta}}$, for which we have

\begin{aligned}\Delta (\mathcal{I}^{-1} \mathbf{p}) = \begin{bmatrix} {\begin{bmatrix} {\begin{bmatrix}\mathcal{I}^{-1} \delta_{rc}\end{bmatrix}}_{r}\end{bmatrix}}_{c} & {\begin{bmatrix}\frac{\partial (\mathcal{I}^{-1})}{\partial \theta_c} \mathbf{p} \\ \end{bmatrix}}_{c}\end{bmatrix} \\ \end{aligned}

The $\delta$ over all rows $r$ and columns $c$ is the identity matrix and we are left with

\begin{aligned}\Delta (\mathcal{I}^{-1} \mathbf{p}) & =\begin{bmatrix}\mathcal{I}^{-1} & {\begin{bmatrix}\frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \mathbf{p} \\ \end{bmatrix}}_{c}\end{bmatrix} \end{aligned} \quad\quad\quad(221)

Next, consider just the $P_\theta$ dependence in the elements of the row vector

\begin{aligned}{\begin{bmatrix}\frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \end{bmatrix}}_r\end{aligned} \quad\quad\quad(222)

We can take derivatives of this, and exploiting the fact that these elements are scalars, so they equal their transpose. Also noting that ${A^{-1}}^\text{T} = {A^\text{T}}^{-1}$, and $\mathcal{I} = \mathcal{I}^\text{T}$, we have

\begin{aligned}\frac{\partial {}}{\partial {P_{\theta_c}}}\left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right) & =\frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} {\begin{bmatrix}\delta_{rc}\end{bmatrix}}_r+\frac{1}{{2}} \left({\begin{bmatrix}\delta_{rc}\end{bmatrix}}_r \right)^\text{T}\mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \\ & =\mathbf{p}^\text{T} \mathcal{I}^{-1} \left( \frac{\partial {}}{\partial {\theta_r}}\frac{1}{{2}} \left( Q + Q^T \right) \right)\mathcal{I}^{-1} {\begin{bmatrix}\delta_{rc}\end{bmatrix}}_r \\ & =\mathbf{p}^\text{T} \mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_r}}\mathcal{I}^{-1} {\begin{bmatrix}\delta_{rc}\end{bmatrix}}_r \\ \end{aligned}

Since we also have $B' B^{-1} + B (B^{-1})' = 0$, for invertable matrixes $B$, this reduces to

\begin{aligned}\frac{\partial {}}{\partial {P_{\theta_c}}}\left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right) & =-\mathbf{p}^\text{T} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_r}}{\begin{bmatrix}\delta_{rc}\end{bmatrix}}_r \\ \end{aligned}

Forming the matrix over all rows $r$, and columns $c$, we get a trailing identity multiplying from the right, and are left with

\begin{aligned}{\begin{bmatrix}\frac{\partial {}}{\partial {P_{\theta_c}}}\left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right)\end{bmatrix}}_{r,c} & ={\begin{bmatrix}-\mathbf{p}^\text{T} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_r}}\end{bmatrix}}_{r} ={\begin{bmatrix}-\frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \mathbf{p}\end{bmatrix}}_{c} \end{aligned} \quad\quad\quad(223)

Okay, getting closer. The only thing left is to consider the remaining $\theta$ dependence of 222, and now want the theta partials of the scalar matrix elements

\begin{aligned}\frac{\partial {}}{\partial {\theta_c}} & \left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right) \\ & =\mathbf{p}^\text{T} \left( \frac{\partial {}}{\partial {\theta_c}}\left( \frac{1}{{2}} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \right)\right)\mathbf{p} \\ & =\mathbf{p}^\text{T} \frac{1}{{2}} \mathcal{I}^{-1} \frac{\partial^2 Q^\text{T}}{\partial \theta_c \partial \theta_r} \mathcal{I}^{-1} \mathbf{p} +\mathbf{p}^\text{T} \frac{1}{{2}} \left( \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} +\mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \frac{\partial {(\mathcal{I}^{-1} )}}{\partial {\theta_c}}\right)\mathbf{p} \\ & =\mathbf{p}^\text{T} \frac{1}{{2}} \mathcal{I}^{-1} \frac{\partial^2 Q^\text{T}}{\partial \theta_c \partial \theta_r} \mathcal{I}^{-1} \mathbf{p} +\mathbf{p}^\text{T} \frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \frac{\partial {\mathcal{I}}}{\partial {\theta_r}} \mathcal{I}^{-1} \mathbf{p} \\ \end{aligned}

There is a slight asymmetry between the first and last terms here that can possibly be eliminated. Using ${B^{-1}}' = -B^{-1} B' B^{-1}$, we can factor out the $\mathcal{I}^{-1}\mathbf{p} = \dot{\boldsymbol{\Theta}}$ terms

\begin{aligned}\frac{\partial {}}{\partial {\theta_c}}\left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right) & =\dot{\boldsymbol{\Theta}}^\text{T} \left(\frac{1}{{2}} \frac{\partial^2 Q^\text{T}}{\partial \theta_c \partial \theta_r} -\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1}\frac{\partial {\mathcal{I}}}{\partial {\theta_r}} \right)\dot{\boldsymbol{\Theta}}\end{aligned}

Is this any better? Maybe a bit. Since we are forming the matrix over all $r,c$ indexes and can assume mixed partial commutation the transpose can be dropped leaving us with

\begin{aligned}{\begin{bmatrix}\frac{\partial {}}{\partial {\theta_c}}\left( \frac{1}{{2}} \mathbf{p}^\text{T} \mathcal{I}^{-1} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \mathcal{I}^{-1} \mathbf{p} \right) \end{bmatrix}}_{r,c} & ={\begin{bmatrix}\dot{\boldsymbol{\Theta}}^\text{T} \left(\frac{1}{{2}} \frac{\partial^2 Q}{\partial \theta_c \partial \theta_r} -\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1}\frac{\partial {\mathcal{I}}}{\partial {\theta_r}} \right)\dot{\boldsymbol{\Theta}}\end{bmatrix}}_{r,c} \end{aligned} \quad\quad\quad(224)

We can now assemble all these individual derivatives

\begin{aligned}\dot{\mathbf{z}} & \approx {{\begin{bmatrix}{\begin{bmatrix}\frac{1}{{2}} \dot{\boldsymbol{\Theta}} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_r- \nabla \phi \\ \dot{\boldsymbol{\Theta}}\end{bmatrix}}}_{{t=0}} +{{\begin{bmatrix}-{\begin{bmatrix}\frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \mathbf{p}\end{bmatrix}}_{c} & {\begin{bmatrix}\dot{\boldsymbol{\Theta}}^\text{T} \left(\frac{1}{{2}} \frac{\partial^2 Q}{\partial \theta_c \partial \theta_r} -\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1}\frac{\partial {\mathcal{I}}}{\partial {\theta_r}} \right)\dot{\boldsymbol{\Theta}}-\frac{\partial^2 \phi}{ \partial \theta_r \partial \theta_c }\end{bmatrix}}_{r,c} \\ \\ \mathcal{I}^{-1} & {\begin{bmatrix}\frac{\partial {(\mathcal{I}^{-1})}}{\partial {\theta_c}} \mathbf{p} \\ \end{bmatrix}}_{c}\end{bmatrix} }}_{{t=0}} \mathbf{z}\end{aligned} \quad\quad\quad(225)

We have both ${\partial {(\mathcal{I}^{-1})}}/{\partial {\theta_k}}$ and ${\partial {\mathcal{I}}}/{\partial {\theta_k}}$ derivatives above, which will complicate things when trying to evaluate this for any specific system. A final elimination of the derivatives of the inverse inertial matrix leaves us with

\begin{aligned}\dot{\mathbf{z}} & \approx {{\begin{bmatrix}{\begin{bmatrix}\frac{1}{{2}} \dot{\boldsymbol{\Theta}} \left(\frac{\partial {Q}}{\partial {\theta_r}}\right)^\text{T} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_r- \nabla \phi \\ \dot{\boldsymbol{\Theta}}\end{bmatrix}}}_{{t=0}} +{{\begin{bmatrix}{\begin{bmatrix}\mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_{c} & {\begin{bmatrix}\dot{\boldsymbol{\Theta}}^\text{T} \left(\frac{1}{{2}} \frac{\partial^2 Q}{\partial \theta_c \partial \theta_r} -\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1}\frac{\partial {\mathcal{I}}}{\partial {\theta_r}} \right)\dot{\boldsymbol{\Theta}}-\frac{\partial^2 \phi}{ \partial \theta_r \partial \theta_c }\end{bmatrix}}_{r,c} \\ \\ \mathcal{I}^{-1} & -{\begin{bmatrix}\mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_{c}\end{bmatrix} }}_{{t=0}} \mathbf{z}\end{aligned} \quad\quad\quad(226)

Single pendulum verification.

Having accumultated this unholy mess of abstraction, lets verify this first against the previous result obtained for the single planar pendulum. Then if that checks out, calculate these matrixes explicitly for the double and multiple pendulum cases. For the single mass pendulum we have

\begin{aligned}Q & = \mathcal{I} = m l^2 \\ \phi & = - m g l \cos\theta\end{aligned} \quad\quad\quad(227)

So all the $\theta$ partials except that of the potential are zero. For the potential we have

\begin{aligned}{\left.{{- \frac{\partial^2 \phi}{\partial^2 \theta}}}\right\vert}_{{0}} = - m g l \cos\theta_0\end{aligned} \quad\quad\quad(229)

\begin{aligned}{\left.{{-\nabla \phi}}\right\vert}_{{0}} = \begin{bmatrix} - m g l \sin\theta_0\end{bmatrix} \end{aligned} \quad\quad\quad(230)

Putting these all together in this simplest application of 226 we have for the linear approximation of a single point mass pendulum about some point in phase space at time zero:

\begin{aligned}\dot{\mathbf{z}} \approx \begin{bmatrix} - m g l \sin\theta_0 \\ \dot{\theta}_0\end{bmatrix} +\begin{bmatrix} 0 & - m g l \cos\theta_0 \\ \frac{1}{{m l^2}} & 0\end{bmatrix} \mathbf{z}\end{aligned} \quad\quad\quad(231)

Excellent. Have not gotten into too much trouble with the math so far. This is consistent with the previous results obtained considering the simple pendulum directly (it actually pointed out an error in the earlier pendulum treatment which is now fixed (I’d dropped the $\dot{\theta}_0$ term)).

Double pendulum explicitly.

For the double pendulum, with $\delta = \theta_1 - \theta_2$, and $M = m_1 + m_2$, we have

\begin{aligned}Q = \begin{bmatrix}M {l_1}^2 & m_2 l_2 l_1 e^{i(\theta_2 - \theta_1)} \\ m_2 l_1 l_2 e^{i(\theta_1 - \theta_2)} & m_2 {l_2}^2 \\ \end{bmatrix}=\begin{bmatrix}M {l_1}^2 & m_2 l_2 l_1 e^{-i\delta} \\ m_2 l_1 l_2 e^{i\delta} & m_2 {l_2}^2 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(232)

\begin{aligned}\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\text{T} \left(\frac{\partial {Q}}{\partial {\theta_1}} \right)^\text{T} \dot{\boldsymbol{\Theta}} & = \frac{1}{{2}} m_2 l_1 l_2 i\dot{\boldsymbol{\Theta}}^\text{T} {\begin{bmatrix}0 & - e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix}}^\text{T}\dot{\boldsymbol{\Theta}} \\ & = \frac{1}{{2}} m_2 l_1 l_2 i\dot{\boldsymbol{\Theta}}^\text{T} \begin{bmatrix}e^{i\delta} \dot{\theta}_2 \\ -e^{-i\delta} \dot{\theta}_1 \end{bmatrix} \\ & = \frac{1}{{2}} m_2 l_1 l_2 i \dot{\theta}_1 \dot{\theta}_2 (e^{i\delta} -e^{-i\delta}) \\ & = -m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin\delta\end{aligned}

The $\theta_2$ derivative is the same but inverted in sign, so we have most of the constant term calculated. We need the potential gradient to complete. Our potential was

\begin{aligned}\phi = - M l_1 g \cos{\theta_1} - m_2 l_2 g \cos{\theta_2}\end{aligned} \quad\quad\quad(233)

\begin{aligned}\nabla \phi =\begin{bmatrix}M l_1 g \sin{\theta_1} \\ m_2 l_2 g \sin{\theta_2}\end{bmatrix}\end{aligned} \quad\quad\quad(234)

Putting things back together we have for the linear approximation of the two pendulum system

\begin{aligned}\dot{\mathbf{z}} = {{\begin{bmatrix}m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) \begin{bmatrix}-1 \\ 1\end{bmatrix} & - g\begin{bmatrix}M l_1 \sin{\theta_1} \\ m_2 l_2 \sin{\theta_2} \\ \end{bmatrix} \\ \dot{\theta}_1 \\ \dot{\theta}_2 \\ \end{bmatrix}}}_{{t=0}}+ A \mathbf{z}\end{aligned} \quad\quad\quad(235)

Where $A$ is still to be determined (from 226).

One of the elements of $A$ are the matrix of potential derivatives. These are

\begin{aligned}\begin{bmatrix}\frac{\partial {\nabla \phi}}{\partial {\theta_1}} & \frac{\partial {\nabla \phi}}{\partial {\theta_2}}\end{bmatrix}=\begin{bmatrix}M l_1 g \cos\theta_1 & 0 \\ 0 & m_2 l_2 g \cos\theta_2\end{bmatrix}\end{aligned} \quad\quad\quad(236)

We also need the inertial matrix and its inverse. These are

\begin{aligned}\mathcal{I}=\begin{bmatrix}M {l_1}^2 & m_2 l_2 l_1 \cos\delta \\ m_2 l_1 l_2 \cos\delta & m_2 {l_2}^2 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(237)

\begin{aligned}\mathcal{I}^{-1}=\frac{1}{{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}}\begin{bmatrix}m_2 {l_2}^2 & -m_2 l_2 l_1 \cos\delta \\ -m_2 l_1 l_2 \cos\delta & M {l_1}^2 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(238)

Since

\begin{aligned}\frac{\partial {Q}}{\partial {\theta_1}} & = m_2 l_1 l_2 i\begin{bmatrix}0 & - e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix} \end{aligned} \quad\quad\quad(239)

We have

\begin{aligned}\frac{\partial {}}{\partial {\theta_1}} \frac{\partial {Q}}{\partial {\theta_1}} & =-m_2 l_1 l_2 \begin{bmatrix}0 & e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix} \\ \frac{\partial {}}{\partial {\theta_2}} \frac{\partial {Q}}{\partial {\theta_1}} & =m_2 l_1 l_2 \begin{bmatrix}0 & e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix} \\ \frac{\partial {}}{\partial {\theta_1}} \frac{\partial {Q}}{\partial {\theta_2}} & =m_2 l_1 l_2 \begin{bmatrix}0 & e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix} \\ \frac{\partial {}}{\partial {\theta_2}} \frac{\partial {Q}}{\partial {\theta_2}} & =-m_2 l_1 l_2 \begin{bmatrix}0 & e^{-i\delta} \\ e^{i\delta} & 0\end{bmatrix} \end{aligned} \quad\quad\quad(240)

and the matrix of derivatives becomes

\begin{aligned}\frac{1}{{2}}\dot{\boldsymbol{\Theta}}^\text{T} \frac{\partial {}}{\partial {\theta_c}}\frac{\partial {Q}}{\partial {\theta_r}} \dot{\boldsymbol{\Theta}}=m_2 l_1 l_2 \dot{\theta_1} \dot{\theta_2} \cos(\theta_1 - \theta_2)\begin{bmatrix}-1 & 1 \\ 1 & -1\end{bmatrix} \end{aligned} \quad\quad\quad(244)

For the remaining two types of terms in the matrix $A$ we need $\mathcal{I}^{-1} {\partial {\mathcal{I}}}/{\partial {\theta_k}}$. The derivative of the inertial matrix is

\begin{aligned}\frac{\partial {\mathcal{I}}}{\partial {\theta_k}}=-m_2 l_1 l_2 (\delta_{k1} - \delta_{k2})\begin{bmatrix}0 & \sin\delta \\ \sin\delta & 0\end{bmatrix}\end{aligned} \quad\quad\quad(245)

Computing the product

\begin{aligned}\mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_k}} & =\frac{ -m_2 l_1 l_2 (\delta_{k1} - \delta_{k2}) }{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}m_2 {l_2}^2 & -m_2 l_2 l_1 \cos\delta \\ -m_2 l_1 l_2 \cos\delta & M {l_1}^2 \\ \end{bmatrix}\begin{bmatrix}0 & \sin\delta \\ \sin\delta & 0\end{bmatrix} \\ & =\frac{ -m_2 l_1 l_2 (\delta_{k1} - \delta_{k2}) \sin\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}-m_2 l_2 l_1 \cos\delta & m_2 {l_2}^2 \\ M {l_1}^2 & -m_2 l_1 l_2 \cos\delta \\ \end{bmatrix}\end{aligned}

We want the matrix of $\mathcal{I}^{-1} {\partial {\mathcal{I}}}/{\partial {\theta_c}} \dot{\boldsymbol{\Theta}}$ over columns $c$, and this is

\begin{aligned}{\begin{bmatrix}\mathcal{I}^{-1} {\partial {\mathcal{I}}}/{\partial {\theta_c}} \dot{\boldsymbol{\Theta}}\end{bmatrix}}_{c}=\frac{ m_2 l_1 l_2 \sin\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}m_2 l_2 l_1 \cos\delta \dot{\theta}_1 - m_2 {l_2}^2 \dot{\theta}_2 & -m_2 l_2 l_1 \cos\delta \dot{\theta}_1 + m_2 {l_2}^2 \dot{\theta}_2 \\ -M {l_1}^2 \dot{\theta}_1 +m_2 l_1 l_2 \cos\delta \dot{\theta}_2 & M {l_1}^2 \dot{\theta}_1 -m_2 l_1 l_2 \cos\delta \dot{\theta}_2\end{bmatrix}\end{aligned} \quad\quad\quad(246)

Very messy. Perhaps it would be better not even bothering to expand this explicitly? The last term in the matrix $A$ is probably no better. For that we want

\begin{aligned}-\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_r}} & =\frac{ -{m_2}^2 {l_1}^2 {l_2}^2 (\delta_{c1} - \delta_{c2}) (\delta_{r1} - \delta_{r2}) \sin^2\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}\begin{bmatrix}-m_2 l_2 l_1 \cos\delta & m_2 {l_2}^2 \\ M {l_1}^2 & -m_2 l_1 l_2 \cos\delta \\ \end{bmatrix} \\ & =\frac{ -{m_2}^2 {l_1}^2 {l_2}^2 (\delta_{c1} - \delta_{c2}) (\delta_{r1} - \delta_{r2}) \sin^2\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}M {l_1}^2 & -m_2 l_1 l_2 \cos\delta \\ -m_2 l_2 l_1 \cos\delta & m_2 {l_2}^2 \\ \end{bmatrix} \\ \end{aligned}

With a sandwich of this between $\dot{\boldsymbol{\Theta}}^\text{T}$ and $\dot{\boldsymbol{\Theta}}$ we are almost there

\begin{aligned}-\dot{\boldsymbol{\Theta}}^\text{T}\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_r}}\dot{\boldsymbol{\Theta}} & =\frac{ -{m_2}^2 {l_1}^2 {l_2}^2 (\delta_{c1} - \delta_{c2}) (\delta_{r1} - \delta_{r2}) \sin^2\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\left( M {l_1}^2 {\dot{\theta}_1}^2 -2m_2 l_1 l_2 \cos\delta \dot{\theta_1} \dot{\theta}_2 ++ m_2 {l_2}^2 {\dot{\theta}_2}^2 \right)\end{aligned}

we have a matrix of these scalars over $r,c$, and that is

\begin{aligned}{\begin{bmatrix}-\dot{\boldsymbol{\Theta}}^\text{T}\frac{\partial {\mathcal{I}}}{\partial {\theta_c}} \mathcal{I}^{-1} \frac{\partial {\mathcal{I}}}{\partial {\theta_r}}\dot{\boldsymbol{\Theta}}\end{bmatrix}}_{rc} & =\frac{ {m_2}^2 {l_1}^2 {l_2}^2 \sin^2\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\left(M {l_1}^2 {\dot{\theta}_1}^2 - 2 m_2 l_1 l_2 \cos\delta \dot{\theta_1} \dot{\theta}_2 + m_2 {l_2}^2 {\dot{\theta}_2}^2 \right)\begin{bmatrix}-1 & 1 \\ 1 & -1\end{bmatrix} \end{aligned} \quad\quad\quad(247)

Putting all the results for the matrix $A$ together is going to make a disgusting mess, so lets summarize in block matrix form

\begin{aligned}A & = {{\begin{bmatrix}B & C \\ \mathcal{I}^{-1} & -B\end{bmatrix}}}_{{t=0}} \\ B & =\frac{ m_2 l_1 l_2 \sin\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\begin{bmatrix}m_2 l_2 l_1 \cos\delta \dot{\theta}_1 - m_2 {l_2}^2 \dot{\theta}_2 & -m_2 l_2 l_1 \cos\delta \dot{\theta}_1 + m_2 {l_2}^2 \dot{\theta}_2 \\ -M {l_1}^2 \dot{\theta}_1 +m_2 l_1 l_2 \cos\delta \dot{\theta}_2 & M {l_1}^2 \dot{\theta}_1 -m_2 l_1 l_2 \cos\delta \dot{\theta}_2\end{bmatrix} \\ C & =\left(m_2 l_1 l_2 \dot{\theta_1} \dot{\theta_2} \cos \delta+\frac{ {m_2}^2 {l_1}^2 {l_2}^2 \sin^2\delta}{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}\left(M {l_1}^2 {\dot{\theta}_1}^2 - 2 m_2 l_1 l_2 \cos\delta \dot{\theta_1} \dot{\theta}_2 + m_2 {l_2}^2 {\dot{\theta}_2}^2 \right)\right)\begin{bmatrix}-1 & 1 \\ 1 & -1\end{bmatrix} \\ & \qquad+\begin{bmatrix}M l_1 g \cos\theta_1 & 0 \\ 0 & m_2 l_2 g \cos\theta_2\end{bmatrix} \\ \mathcal{I}^{-1} & =\frac{1}{{{l_1}^2 {l_2}^2 m_2 (M - m_2 \cos^2\delta)}}\begin{bmatrix}m_2 {l_2}^2 & -m_2 l_2 l_1 \cos\delta \\ -m_2 l_1 l_2 \cos\delta & M {l_1}^2 \\ \end{bmatrix} \\ b & =\begin{bmatrix}m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1 - \theta_2) \begin{bmatrix}-1 \\ 1\end{bmatrix} & - g\begin{bmatrix}M l_1 \sin{\theta_1} \\ m_2 l_2 \sin{\theta_2} \\ \end{bmatrix} \\ \dot{\theta}_1 \\ \dot{\theta}_2 \\ \end{bmatrix}\end{aligned}

where these are all related by the first order matrix equation

\begin{aligned}\frac{d\mathbf{z}}{dt} & = {\left.{\mathbf{b}}\right\vert}_{{t=0}} + {\left.{{A}}\right\vert}_{{t=0}} \mathbf{z}\end{aligned} \quad\quad\quad(248)

Wow, even to just write down the equations required to get a linear approximation of the two pendulum system is horrendously messy, and this isn’t even trying to solve it. Numerical and or symbolic computation is really called for here. If one elected to do this numerically, which looks pretty much mandatory since the analytic way didn’t turn out to be simple even for just the two pendulum system, then one is probably better off going all the way back to 215 and just calculating the increment for the tragectory using a very small time increment, and do this repeatedly (i.e. do a zeroth order numerical procedure instead of the first order which turns out much more complicated).

Posted in Math and Physics Learning. | Tagged: , , , | 4 Comments »

scope conflict between member variable and loop initializer.

Posted by peeterjoot on October 19, 2009

Here’s a scoping gotcha that bit me today … scratched my head over this one for a bit in the debugger before I figured it out.

I was working on some code that had some uninitialized variables, one of which was like so:

unsigned i ;
// pages of junk
for ( i = 0 ; i < blah ; i++ )
{
}
...
for ( i = 0, m_confused = 0 ; i < totalIterations ; i++ ) // m_confused is a member variable.
{
}


There were some other uninitialized variables, and I gave these all explicit invalid or zero values to make sure they didn’t have stack garbage in them. For the loop counter i, I just moved it into the for () declaration like so:

for ( unsigned i = 0 ; i < blah ; i++ )
{
}
...
for ( unsigned i = 0, m_confused = 0 ; i < totalIterations ; i++ )
{
if ( stuff )
{
m_confused = N ;
break ;
}
}


Did you spot my mistake? I didn’t when I made the change, and it was so minor that I forgot I’d even made it by the time I tried out the code.

After the loop termination my member variable m_confused was mysteriously set to zero, despite the loop termination condition in the non-error codepath requiring that a value be assigned.

Turns out that the compiler (g++ 4.2.3 in this case) after I “fixed” the code to remove the possibility of referencing the loop counter inappropriately interpreted this as me declaring a loop scoped variable m_confused. That, and not this->m_confused got my update, so once it went out of scope I was left with zero in my class. Unfriendly of the compiler to not give me a warning about this even with -Wall (the older gcc 4.1.2 compiler does).

What are the lessons to be learned?

• don’t fix unbroken code.  This unfortunately takes much more self discipline than I have!
• no more than one statement per line.  Keeping the code simple, should give less chance for you to misinform the compiler what you want to do.  There wasn’t any good reason that the original code needed to initialize both of these in the for loop.  The ‘m_confused = 0’ statement could have just as easily been on it’s own line.  Ironically, this particular value was already zeroed out in other code, so it didn’t even have to be in the for loop explicitly at all.
• always initialize all stack variables no matter what at the point of declaration.  Even the simple change to make it obvious that it is initialized by restricting the usage scope can be screwed up.
• have good code coverage testing.  If I had been unlucky enough to have made this “maintainability-fix” in a non-mainline codepath, I’d have made a very unfriendly alteration for somebody else (or myself) to later debug.

Ontario Science Center encouraging teen mothers?

Posted by peeterjoot on October 12, 2009

Doesn’t this picture speak for itself?

Are you or your child between the ages of 7 - 17

Hamiltonian treatment of rigid Spherical pendulum.

Posted by peeterjoot on October 11, 2009

For the spherical rigid pendulum of length $l$, we have for the distance above the lowest point

\begin{aligned}h = l (1 + \cos\theta)\end{aligned} \quad\quad\quad(157)

(measuring $\theta$ down from the North pole as conventional). The Lagrangian is therefore

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m l^2 (\dot{\theta}^2 + \sin^2\theta \dot{\phi}^2) - m g l (1 + \cos\theta)\end{aligned} \quad\quad\quad(158)

We can drop the constant term, using the simpler Lagrangian

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m l^2 (\dot{\theta}^2 + \sin^2\theta \dot{\phi}^2) - m g l \cos\theta\end{aligned} \quad\quad\quad(159)

To express the Hamiltonian we need first the conjugate momenta, which are

\begin{aligned}P_\theta &= \frac{\partial {\mathcal{L}}}{\partial {\dot{\theta}}} = m l^2 \dot{\theta} \\ P_\phi &= \frac{\partial {\mathcal{L}}}{\partial {\dot{\phi}}} = m l^2 \sin^2\theta \dot{\phi}\end{aligned} \quad\quad\quad(160)

We can now write the Hamiltonian

\begin{aligned}H = \frac{1}{{2 m l^2}} \left({P_\theta}^2 + \frac{1}{{\sin^2\theta}} {P_\phi}^2\right) + m g l \cos\theta\end{aligned} \quad\quad\quad(162)

Before going further one sees that there is going to be trouble where $\sin\theta = 0$. Curiously, this is at the poles, the most dangling position and the upright. The south pole is the usual point where we solve the planar pendulum problem using the harmonic oscillator approximation, so it is somewhat curious that the energy of the system appears to go undefined at this point where the position is becoming more defined. It seems almost like a quantum uncertainty phenomena until one realizes that the momentum conjugate to $\phi$ is itself proportional to $\sin^2 \theta$. By expressing the energy in terms of this $P_\phi$ momentum we have to avoid looking at the poles for a solution to the equations. If we go back to the Lagrangian and the Euler-Lagrange equations, this point becomes perfectly tractable since we are no longer dividing through by $\sin^2\theta$.

Examining the polar solutions is something to return to. For now, let’s avoid that region. For regions where $\sin\theta$ is nicely non-zero, we get for the Hamiltonian equations

\begin{aligned}\frac{\partial {H}}{\partial {P_\phi}} &= \dot{\phi} = \frac{1}{{ m l^2 \sin^2 \theta}} P_\phi \\ \frac{\partial {H}}{\partial {P_\theta}} &= \dot{\theta} = \frac{1}{{ m l^2 }} P_\theta \\ \frac{\partial {H}}{\partial {\phi}} &= -\dot{P}_\phi = 0 \\ \frac{\partial {H}}{\partial {\theta}} &= -\dot{P}_\theta = -\frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 - m g l \sin\theta\end{aligned} \quad\quad\quad(163)

These now expressing the dynamics of the system. The first two equations are just the definitions of the canonical momenta that we started with using the Lagrangian. Not surprisingly, but unfortunate, we have a non-linear system here like the planar rigid pendulum, so despite this being one of the most simple systems it does not look terribly tractable. What would it take to linearize this system of equations?

Lets write the state space vector for the system as

\begin{aligned}\mathbf{x} =\begin{bmatrix}P_\theta \\ \theta \\ P_\phi \\ \phi \\ \end{bmatrix}\end{aligned} \quad\quad\quad(167)

lets also suppose that we are interested in the change to the state vector in the neighborhood of an initial state

\begin{aligned}\mathbf{x} =\begin{bmatrix}P_\theta \\ \theta \\ P_\phi \\ \phi \\ \end{bmatrix}={\begin{bmatrix}P_\theta \\ \theta \\ P_\phi \\ \phi \\ \end{bmatrix}}_0+ \mathbf{z}\end{aligned} \quad\quad\quad(168)

The Hamiltonian equations can then be written

\begin{aligned}\frac{d\mathbf{z}}{dt} &=\begin{bmatrix}\frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 + m g l \sin\theta \\ \frac{1}{{ m l^2 }} P_\theta \\ 0 \\ \frac{1}{{ m l^2 \sin^2 \theta}} P_\phi \\ \end{bmatrix}\end{aligned} \quad\quad\quad(169)

Getting away from the specifics of this particular system is temporarily helpful. We have a set of equations that we wish to calculate a linear approximation for

\begin{aligned}\frac{dz_\mu}{dt} = A_\mu(x_\nu) \approx A_\mu(\mathbf{x}_0) + \sum_\alpha {\left. \frac{\partial {A_\mu}}{\partial {x_\alpha}} \right\vert}_{\mathbf{x}_0} z_\alpha\end{aligned} \quad\quad\quad(170)

Our linear approximation is thus

\begin{aligned}\frac{d\mathbf{z}}{dt} &\approx{\begin{bmatrix}\frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 + m g l \sin\theta \\ \frac{1}{{ m l^2 }} P_\theta \\ 0 \\ \frac{1}{{ m l^2 \sin^2 \theta}} P_\phi \\ \end{bmatrix}}_0+{\begin{bmatrix}0 & -\frac{{P_\phi}^2 (1 + 2 \cos^2 \theta)}{m l^2 \sin^4 \theta} +m g l \cos\theta & \frac{2 \cos\theta}{ m l^2 \sin^3 \theta} P_\phi & 0 \\ \frac{1}{{m l^2}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & \frac{-2 P_\phi \cos\theta}{m l^2 \sin^3\theta} & \frac{1}{{ m l^2 \sin^2 \theta}} & 0 \\ \end{bmatrix}}_0\mathbf{z}\end{aligned} \quad\quad\quad(171)

Now, this is what we get blindly trying to set up the linear approximation of the state space differential equation. We see that the cyclic coordinate $\phi$ leads to a bit of trouble since no explicit $\phi$ dependence in the Hamiltonian makes the resulting matrix factor non-invertible. It appears that we would be better explicitly utilizing this cyclic coordinate to note that $P_\phi = \text{constant}$, and to omit this completely from the state vector. Our equations in raw form are now

\begin{aligned}\dot{\theta} &= \frac{1}{{ m l^2 }} P_\theta \\ \dot{P}_\theta &= \frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 + m g l \sin\theta \\ \dot{\phi} &= \frac{1}{{ m l^2 \sin^2 \theta}} P_\phi \end{aligned} \quad\quad\quad(172)

We can treat the $\phi$ dependence later once we have solved for $\theta$. That equation to later solve is just this last

\begin{aligned}\dot{\phi} &= \frac{1}{{ m l^2 \sin^2 \theta}} P_\phi \end{aligned} \quad\quad\quad(175)

This integrates directly, presuming $\theta = \theta(t)$ is known, and we have

\begin{aligned}\phi - \phi(0) &= \frac{P_\phi}{ m l^2} \int_0^t \frac{d\tau}{\sin^2 \theta(\tau)} \end{aligned} \quad\quad\quad(176)

Now the state vector and its perturbation can be redefined omitting all but the $\theta$ dependence. Namely

\begin{aligned}\mathbf{x} =\begin{bmatrix}P_\theta \\ \theta \\ \end{bmatrix}\end{aligned} \quad\quad\quad(177)

\begin{aligned}\mathbf{x} =\begin{bmatrix}P_\theta \\ \theta \\ \end{bmatrix}={\begin{bmatrix}P_\theta \\ \theta \\ \end{bmatrix}}_0+ \mathbf{z}\end{aligned} \quad\quad\quad(178)

We can now write the remainder of this non-linear system as

\begin{aligned}\frac{d\mathbf{z}}{dt} &=\begin{bmatrix}\frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 + m g l \sin\theta \\ \frac{1}{{ m l^2 }} P_\theta \\ \end{bmatrix}\end{aligned} \quad\quad\quad(179)

and make the linear approximation around $\mathbf{x}_0$ as

\begin{aligned}\frac{d\mathbf{z}}{dt} &\approx{\begin{bmatrix}\frac{\cos\theta}{ m l^2 \sin^3 \theta} {P_\phi}^2 + m g l \sin\theta \\ \frac{1}{{ m l^2 }} P_\theta \\ \end{bmatrix}}_0+{\begin{bmatrix}0 & -\frac{{P_\phi}^2 (1 + 2 \cos^2 \theta)}{m l^2 \sin^4 \theta} +m g l \cos\theta \\ \frac{1}{{m l^2}} & 0 \\ \end{bmatrix}}_0\mathbf{z}\end{aligned} \quad\quad\quad(180)

This now looks a lot more tractable, and is in fact exactly the same form now as the equation for the linearized planar pendulum. The only difference is the normalization required to switch to less messy dimensionless variables. The main effect of allowing the trajectory to have a non-planar component is a change in the angular frequency in the $\theta$ dependent motion. That frequency will no longer be $\sqrt{{\left\lvert{\cos\theta_0}\right\rvert} g/l}$, but also has a $P_\phi$ and other more complex trigonometric $\theta$ dependencies. It also appears that we can probably have hyperbolic or trigonometric solutions in the neighborhood of any point, regardless of whether it is a northern hemispherical point or a southern one. In the planar pendulum the unambiguous sign of the matrix terms led to hyperbolic only above the horizon, and trigonometric only below.

What’s in a name. The meaning of Joot.

Posted by peeterjoot on October 10, 2009

I’ve been listening to Estonian podcasts. The most recent series concerned the guilds, merchants and so on starting in the Middle Ages and going into the 17th Century. I had heard a couple of these, and the subject material was of very little interest to me until I hit a couple of them which sheds some light about the origin of our family name.

I don’t necessarily want to transcribe the episodes, but will give a summary of what I heard.

Basically a “joot” was a festive occasion ( the use of the word in this sense is not current used in the language.) Strictly speaking these “joots” were feasts usually held in guild ( trade or merchant associations )halls, where there was plenty to eat drink. Sometimes dancing was an element in these “joots”.

There were three major events during the year where “ joots ” played an important part in festivities. These were the times around Christmas, Easter, and the coming of summer in the month of May. The Christmas and Easter events were usually two weeks long, and the one in May one week. There were other events during the year which were not as long, but all of them ended with a “joot”.

The two week events involved parades, musicians ( a lot of drummers and (bag-pipers )dancing, singing, presentations of all kinds of events – street dancing and singing, sword, spear and torch dances done by men, costumes ( one favourite one was that of the Devil). These were much in the spirit of a carnival and the festivities were always followed by a “joot”

The May celebrations which usually lasted for a week, where a “prince of spring” was chosen always ended up with a “joot”. Some of these involved burning of effigies etc.

Some of the shorter celebrations or events involved target shooting ( first with cross-bows and later with rifles ), horse races, and contests of various kinds always followed by a “joot” where there was plenty to eat and drink.

1.)Although the word “joot” is not today used in the mediaeval sense of a gala “state” celebration, it still today carries on the connotation of a festive event which consists of a get-together of which drinking is a definite part. I have never heard the word or seen the word in print with this meaning, but it does exist in a BIG dictionary.( The action word of “jooma” means to drink, and a “joot” would describe the event where drinking took place.

2.)There is another word “jootma” which means to solder. A jootraud would be a soldering iron. ( Raud means iron )

3.) There is a combining word form “jootraha” (“raha” means money) which means essentially a tip such that you would give for good service in a restaurant. There are apparently some dictionaries in Estonian which trace back the words to their source, but I don’t own one of them. I would suspect that the derivation of “jootraha” goes back to number 1.). This would probably go to the custom of tipping of what I would assume to be the serving wenches for good service or the hope of good future service. Alternately and more probably, the giving of tips had its origin in the participants of “joots” giving money to the organizers of these in the way of fixed “donations” or donations in accordance to their ability to pay. So its joining in to pay for the celebrations.

This covers the variations to some degree in the name.

Line element for Spherical Polar unit vectors using Geometric Algebra.

Posted by peeterjoot on October 9, 2009

The line element for the particle moving on a spherical surface can be calculated by calculating the derivative of the spherical polar unit vector $\hat{\mathbf{r}}$

\begin{aligned}\frac{d\hat{\mathbf{r}}}{dt} = \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \frac{d\phi}{dt} +\frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \frac{d\theta}{dt}\end{aligned} \quad\quad\quad(17)

than taking the magnitude of this vector. We can start either in coordinate form

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta \cos\phi + \mathbf{e}_2 \sin\theta \sin\phi\end{aligned} \quad\quad\quad(18)

or, instead do it the fun way, first grouping this into a complex exponential form. Writing $i = \mathbf{e}_1 \mathbf{e}_2$, the first factorization is

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta e^{i\phi} \end{aligned} \quad\quad\quad(19)

The unit vector $\boldsymbol{\rho} = \mathbf{e}_1 e^{i\phi}$ lies in the $x,y$ plane perpendicular to $\mathbf{e}_3$, so we can form the unit bivector $\mathbf{e}_3\boldsymbol{\rho}$ and further factor the unit vector terms into a $\cos + i \sin$ form

\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_3 \cos\theta + \mathbf{e}_1 \sin\theta e^{i\phi} \\ &= \mathbf{e}_3 (\cos\theta + \mathbf{e}_3 \boldsymbol{\rho} \sin\theta) \\ \end{aligned}

This allows the spherical polar unit vector to be expressed in complex exponential form (really a vector-quaternion product)

\begin{aligned}\hat{\mathbf{r}} = \mathbf{e}_3 e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} = e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \mathbf{e}_3\end{aligned} \quad\quad\quad(20)

Now, calcuating the unit vector velocity, we get

\begin{aligned}\frac{d\hat{\mathbf{r}}}{dt} &= \mathbf{e}_3 \mathbf{e}_3 \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \dot{\theta} + \mathbf{e}_1 \mathbf{e}_1 \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \\ &= \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \left(\dot{\theta} + e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho} \sin\theta e^{-i\phi} \mathbf{e}_2 \dot{\phi}\right) \\ &= \left( \dot{\theta} + \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \right) e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho}\end{aligned}

The last two lines above factor the $\boldsymbol{\rho}$ vector and the $e^{\mathbf{e}_3 \boldsymbol{\rho} \theta}$ quaternion to the left and to the right in preparation for squaring this to calculate the magnitude.

\begin{aligned}\left( \frac{d\hat{\mathbf{r}}}{dt} \right)^2&=\left\langle{{ \left( \frac{d\hat{\mathbf{r}}}{dt} \right)^2 }}\right\rangle \\ &=\left\langle{{ \left( \dot{\theta} + \mathbf{e}_2 \sin\theta e^{i\phi} \dot{\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \boldsymbol{\rho} \theta} \right) \left(\dot{\theta} + e^{-\mathbf{e}_3 \boldsymbol{\rho} \theta} \boldsymbol{\rho} \sin\theta e^{-i\phi} \mathbf{e}_2 \dot{\phi}\right) }}\right\rangle \\ &=\dot{\theta}^2 + \sin^2\theta \dot{\phi}^2+ \sin\theta \dot{\phi} \dot{\theta}\left\langle{{ \mathbf{e}_2 e^{i\phi} \boldsymbol{\rho} e^{\mathbf{e}_3 \rho \theta}+e^{-\mathbf{e}_3 \rho \theta} \boldsymbol{\rho} e^{-i\phi} \mathbf{e}_2}}\right\rangle \\ \end{aligned}

This last term ($\in \text{span} \{\rho\mathbf{e}_1, \rho\mathbf{e}_2, \mathbf{e}_1\mathbf{e}_3, \mathbf{e}_2\mathbf{e}_3\}$) has only grade two components, so the scalar part is zero. We are left with the line element

\begin{aligned}\left(\frac{d (r\hat{\mathbf{r}})}{dt}\right)^2 = r^2 \left( \dot{\theta}^2 + \sin^2\theta \dot{\phi}^2 \right)\end{aligned} \quad\quad\quad(21)

In retrospect, at least once one sees the answer, it seems obvious. Keeping $\theta$ constant the length increment moving in the plane is $ds = \sin\theta d\phi$, and keeping $\phi$ constant, we have $ds = d\theta$. Since these are perpendicular directions we can add the lengths using the Pythagorean theorem.