[Click here for a PDF of this post with nicer formatting]
(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
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
axis down, and
axis to the left with
we have for the position of the first mass
, at angle
and length 

If the second mass, dangling from this is at an angle
from the
axis, its position is

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

For the second mass

Taking conjugates and multiplying out we have

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

For the second particle, the distance from the horizontal is

So the total distance from the reference point is

We now have the Lagrangian

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

The conjugate momenta that we need for the Hamiltonian are

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

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
, looks like it is something like a moment of inertia tensor. If we call this
, and write

Then the relation between the conjugate momenta in vector form

and the angular velocity vector can be written

Can we write the Lagrangian in terms of
? The first Kinetic term is easy, just

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

Writing
for this
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

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


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

In the expression for
above, it is implied that the matrix is zero for any indexes
, so it would perhaps be better to write explicitly

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

Similarly the
conjugate momentum is

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

Observing that
, and thus
, we now have everything required to express the Hamiltonian in terms of the conjugate momenta

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

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

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

The remaining Hamiltonian equations are
, and the tricky part of evaluating this is going to all reside in the Kinetic term. Diving right in this is

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

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

This gives us for the Hamiltonian equation

If we introduce a phase space position gradients

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

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

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
, and
derivatives of the top
rows of this matrix.
If we write

and the Hamiltonian equations 215 as

Then the linearization, without simplifying or making explicit yet is

For brevity the constant term evaluated at
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
for the matrix of first order derivatives in our Taylor expansion, as in

First, lets do the potential gradient.

Next in terms of complexity is the first order term of
, for which we have

The
over all rows
and columns
is the identity matrix and we are left with

Next, consider just the
dependence in the elements of the row vector

We can take derivatives of this, and exploiting the fact that these elements are scalars, so they equal their transpose. Also noting that
, and
, we have

Since we also have
, for invertable matrixes
, this reduces to

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

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

There is a slight asymmetry between the first and last terms here that can possibly be eliminated. Using
, we can factor out the
terms

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

We can now assemble all these individual derivatives

We have both
and
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

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

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

and for the angular gradient

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:

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
term)).
Double pendulum explicitly.
For the double pendulum, with
, and
, we have


The
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

So, the gradient is

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

Where
is still to be determined (from 226).
One of the elements of
are the matrix of potential derivatives. These are

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


Since

We have

and the matrix of derivatives becomes

For the remaining two types of terms in the matrix
we need
. The derivative of the inertial matrix is

Computing the product

We want the matrix of
over columns
, and this is

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

With a sandwich of this between
and
we are almost there

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

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

where these are all related by the first order matrix equation

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).