Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Spherical polar pendulum for one and multiple masses, and multivector Euler-Lagrange formulation.

Posted by peeterjoot on November 4, 2009

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

Motivation

The planar multiple pendulum problem proved somewhat tractable in the Hamiltonian formulation. Generalizing this to allow for three dimensional motion is a logical next step. Here this is attempted, using a Geometric Algebra scalar plus bivector parametrization of the spherical position of each dangling mass relative to the position of what it is attached to, as in

\begin{aligned}\mathbf{z} = l \mathbf{e}_3 e^{j\theta}\end{aligned} \quad\quad\quad(1)

The exponential is essentially a unit quaternion, rotating the vector l \mathbf{e}_3 from the polar axis to its \theta,\phi angle dependent position. Two sided rotation operators are avoided here by requiring of the unit bivector j = \mathbf{e}_3 \wedge \mathbf{m}, where \mathbf{m} is a vector in the plane of rotation passing through the great circle from \mathbf{e}_3 through \mathbf{z}. Note that we are free to pick \mathbf{m} = \mathbf{e}_1 e^{\mathbf{e}_1 \mathbf{e}_2 \phi}, the unit vector in the x,y plane at angle \phi from the x-axis. When that is done j = \mathbf{e}_3 \mathbf{m} since these are perpendicular. Setting up the Lagrangian in terms of the bivector j instead of the scalar angle \phi will be attempted, since this is expected to have some elegance and will be a fun way to try the problem. This should also provide a concrete example of a multivector Lagrangian in a context much simpler than electromagnetic fields or quantum mechanics.

Note finally that a number of simplifying assumptions will be made. These include use of point masses, zero friction at the pivots and rigid nonspringy massless connecting rods between the masses.

Kinetic energy for the single pendulum case.

Let’s compute derivatives of the unit vector

\begin{aligned}\hat{\mathbf{z}} = \mathbf{e}_3 e^{j\theta} = e^{-j\theta} \mathbf{e}_3\end{aligned} \quad\quad\quad(2)

This can be done with both the left and right factorization of \mathbf{e}_3, and are respectively

\begin{aligned}\dot{\hat{\mathbf{z}}} &= \mathbf{e}_3 \left( j \dot{\theta} e^{j\theta} + \frac{dj}{dt} \sin\theta \right) \\ &=\mathbf{e}_3\begin{bmatrix}j e^{j\theta} & \sin\theta\end{bmatrix}\frac{d}{dt} \begin{bmatrix}\theta \\ j\end{bmatrix}\end{aligned} \quad\quad\quad(3)

\begin{aligned}\dot{\hat{\mathbf{z}}} &= \left( -j \dot{\theta} e^{-j\theta} - \frac{dj}{dt} \sin\theta \right) \mathbf{e}_3 \\ &=\left(\frac{d}{dt} \begin{bmatrix}\theta & -j\end{bmatrix}\right)\begin{bmatrix}-j e^{-j\theta} \\ \sin\theta\end{bmatrix}\mathbf{e}_3\end{aligned} \quad\quad\quad(5)

These derivatives have been grouped into a matrix factors that allow a natural seeming conjugate operation to be defined. That is for a matrix of multivector elements a_{ij}

\begin{aligned}A =\begin{bmatrix}a_{ij}\end{bmatrix}\end{aligned} \quad\quad\quad(7)

define a conjugate matrix, as the transpose of the reversed elements

\begin{aligned}A^\dagger \equiv\begin{bmatrix}\tilde{a}_{ji}\end{bmatrix}\end{aligned} \quad\quad\quad(8)

With this definition, plus two helpers

\begin{aligned}\boldsymbol{\Theta} &\equiv\begin{bmatrix}\theta \\  j\end{bmatrix} \\ R &= \begin{bmatrix}j e^{j\theta} & \sin\theta\end{bmatrix}\end{aligned} \quad\quad\quad(9)

Our velocity becomes

\begin{aligned}{\dot{\hat{\mathbf{z}}}}^2 &= {\dot{\boldsymbol{\Theta}}}^\dagger R^\dagger R \dot{\boldsymbol{\Theta}} \\ \end{aligned} \quad\quad\quad(11)

Explicitly, expanding the inner matrix product we can write

\begin{aligned}Q &\equiv R^\dagger R \\ &=\begin{bmatrix}1 & -j e^{-j\theta} \sin\theta \\ j e^{j\theta} \sin\theta & \sin^2 \theta\end{bmatrix} \\ {\dot{\hat{\mathbf{z}}}}^2 &= {\dot{\boldsymbol{\Theta}}}^\dagger Q \dot{\boldsymbol{\Theta}} \\ \end{aligned} \quad\quad\quad(13)

This is a slightly unholy mix of geometric and matrix algebra, but it works to compactly express the velocity dependence. Observe that this inner matrix Q = Q^\dagger, so it is Hermitian with this definition of conjugation.

Our Lagrangian for the one particle pendulum, measuring potential energy from the horizontal, is then

\begin{aligned}\mathcal{L} = {\dot{\boldsymbol{\Theta}}}^\dagger \frac{1}{{2}} m l^2 Q \dot{\boldsymbol{\Theta}} - m g l \cos\theta\end{aligned} \quad\quad\quad(17)

We also have a mechanism that should generalize fairly easily to the two or many pendulum cases too.

Before continuing, it should be noted that there were assumptions made in this energy expression derivation that are not reflected in the Lagrangian above. One of these was the unit bivector assumption for j, as well as a \mathbf{e}_3 containment assumption for the plane this represents (\mathbf{e}_3 \wedge j = 0). So for completeness we should probably add to the Lagrangian above some Lagrange multiplier enforced constraints

\begin{aligned}\lambda ( {j}^2 + 1 ) + \alpha \dot{c} ( \mathbf{e}_3 \wedge j )\end{aligned} \quad\quad\quad(18)

Here \alpha has got to be a trivector multiplier for the Lagrangian to be a scalar. Can we get away with omitting these constraints?

Two and multi particle case.

Having constructed a way that can express the velocity of a single spherical pendulum in a tidy way, we can move on to consider the multiple pendulum case as shown in figure (\ref{fig:sPolarMultiPendulum:pendulumDouble})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.4\textheight]{pendulumDouble}
\caption{Double spherical pendulum.}
\end{figure}

There are two bivectors depicted, j_1 and j_2 representing oriented planes passing through great circles from a local polar axis (in direction \mathbf{e}_3). Let the positions of the respective masses be z_1 and z_2, where each mass is connected by a rigid massless rod of length l_1 and l_2 respectively. The masses are rotated by angles \theta_1 and \theta_2 in the planes j_1 and j_2 from an initial direction of \mathbf{e}_3. We can express the position of the second mass as

\begin{aligned}\mathbf{z}_2 = \mathbf{z}_1 + \mathbf{e}_3 e^{j_2 \theta_2}\end{aligned} \quad\quad\quad(19)

We can use the same factorization as previously used for the single mass case and write for our collection of angular velocities

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

Using this the total Kinetic energy is

\begin{aligned}K &= {\dot{\boldsymbol{\Theta}}}^\dagger \frac{1}{{2}} Q \dot{\boldsymbol{\Theta}} \\ R_1 &=\begin{bmatrix}l_1 j_1 e^{j_1\theta_1} & l_1 \sin\theta_1 &0 &0\end{bmatrix} \\ R_2 &=\begin{bmatrix} l_1 j_1 e^{j_1\theta_1} &  l_1 \sin\theta_1 & l_2 j_2 e^{j_2\theta_2} &  l_2 \sin\theta_2 \end{bmatrix} \\ Q &=m_1 {R_1}^\dagger R_1+m_2 {R_2}^\dagger R_2\end{aligned} \quad\quad\quad(21)

Notation has been switched slightly from the single mass case, and the m l^2 factor is now incorporated directly into Q for convenience.

An expansion of Q is essentially one of block matrix multiplication (where we already have to be careful with order of operations as we do for the geometric product elements themselves). We have something like

\begin{aligned}R_1 &= \begin{bmatrix}A_1 & 0\end{bmatrix} \\ {R_1}^\dagger &= \begin{bmatrix}{A_1}^\dagger \\ 0\end{bmatrix} \\ R_2 &= \begin{bmatrix}A_1 & A_2\end{bmatrix} \\ {R_2}^\dagger &= \begin{bmatrix}{A_1}^\dagger \\ {A_2}^\dagger \\ \end{bmatrix} \end{aligned} \quad\quad\quad(25)

We have for the products

\begin{aligned}{R_1}^\dagger R_1 &=\begin{bmatrix}{A_1}^\dagger A_1 & 0 \\ 0 & 0\end{bmatrix} \\ {R_2}^\dagger R_2 &=\begin{bmatrix}{A_1}^\dagger A_1 ^\dagger A_2 \\ {A_2}^\dagger A_1 ^\dagger A_2 \\ \end{bmatrix} \end{aligned} \quad\quad\quad(29)

So our quadratic form matrix is

\begin{aligned}Q =\begin{bmatrix}(m_1 + m_2) {A_1}^\dagger A_1 & m_2 {A_1}^\dagger A_2 \\ m_2 {A_2}^\dagger A_1 & m_2 {A_2}^\dagger A_2 \\ \end{bmatrix} \end{aligned} \quad\quad\quad(31)

In general for the multiple particle case this is

\begin{aligned}Q &={\begin{bmatrix}\left(\sum_{k=\max(r,c)}^N m_k \right){A_r}^\dagger A_c\end{bmatrix}}_{rc} \\ A_k &= l_k \begin{bmatrix}j_k e^{j_k \theta_k} & \sin\theta_k\end{bmatrix}\end{aligned} \quad\quad\quad(32)

Expanded explicitly this is

\begin{aligned}Q &={\begin{bmatrix}\left(\sum_{k=\max(r,c)}^N m_k \right) l_r l_c\begin{bmatrix}- j_r e^{-j_r \theta_r} j_c e^{j_c\theta_c} & - j_r e^{-j_r \theta_r} \sin\theta_c \\ j_c e^{j_c \theta_c} \sin\theta_r & \sin\theta_r \sin\theta_c\end{bmatrix}\end{bmatrix}}_{rc}\end{aligned} \quad\quad\quad(34)

Observe that the order of products in this expansion is specifically ordered, since the j_c and j_r bivectors do not necessarily commute.

The potential in the multiple particle case is also fairly straightforward to compute. Consider the two particle case to illustrate the pattern. Using the lowest point as the potential reference we have

\begin{aligned}\phi' = g \sum m_i h_i= m_1 r_1 (1 + \cos\theta_1) + m_2 \left( r_1(1 + \cos\theta_1) + r_2( 1 + \cos\theta_2) \right)\end{aligned} \quad\quad\quad(35)

Alternately, dropping all the constant terms (using the horizon as the potential reference) we have for the general case

\begin{aligned}\phi = g \sum_i \left( \sum_{k=i}^N m_k \right) r_i \cos\theta_i\end{aligned} \quad\quad\quad(36)

Lets collect all the bits and pieces now for the multiple pendulum Lagrangian now, repeating for coherency, and introducing a tiny bit more notation (mass sums and block angular velocity matrices) for convenience

\begin{aligned}\mathcal{L} &= K - \phi + \sum_k \lambda_k ({j_k}^2 + 1 ) \\ \boldsymbol{\Theta}_i &=\begin{bmatrix}\theta_i \\ j_i\end{bmatrix} \\ \boldsymbol{\Theta} &={\begin{bmatrix}\boldsymbol{\Theta}_r\end{bmatrix}}_{r} \\ \mu_i &=\sum_{k=i}^N m_k \\ Q &={\begin{bmatrix}\mu_{\max(r,c)}l_r l_c\begin{bmatrix}- j_r e^{-j_r \theta_r} j_c e^{j_c\theta_c} & - j_r e^{-j_r \theta_r} \sin\theta_c \\ j_c e^{j_c \theta_c} \sin\theta_r & \sin\theta_r \sin\theta_c\end{bmatrix}\end{bmatrix}}_{rc} \\ K &=\frac{1}{{2}} {\dot{\boldsymbol{\Theta}}}^\dagger Q \dot{\boldsymbol{\Theta}} \\ \phi &=g \sum_{i=1}^N \mu_i r_i \cos\theta_i\end{aligned} \quad\quad\quad(37)

Building up to the Multivector Euler-Lagrange equations.

Rather than diving right into any abstract math, lets consider a few specific examples of multivector Lagrangians to develop some comfort with non-scalar Lagrangian functions. The generalized “coordinates” in the Lagrangian for the spherical pendulum problem being considered include include bivectors, and we don’t know how to evaluate the Euler Lagrange equations for anything but scalar generalized coordinates.

The goal is to develop Euler-Lagrange equations that can handle a Lagrangian for these more general functions, but until we figure out how to do that, we can at least tackle the problem using variation of the action around a stationary solution.

A first example to build intuition.

To help understand what we have to do, lets consider the very simplest bivector parametrized Lagrangian, that of a spherical pendulum constrained (perhaps by a track or a surface) of moving only in a ring. This is shown pictorially in figure (\ref{fig:sPolarMultiPendulum:pendulumPolarCircular})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.4\textheight]{pendulumPolarCircular}
\caption{Circularly constrained spherical pendulum.}
\end{figure}

The potential energy is fixed on this surface, so our Lagrangian is purely kinetic

\begin{aligned}\mathcal{L} = -\frac{1}{{2}} m l^2 \sin^2 \theta_0 {j'}^2\end{aligned} \quad\quad\quad(44)

We’d like to directly vary the action for the Lagrangian around a stationary solution

\begin{aligned}S = \int \mathcal{L} dt\end{aligned} \quad\quad\quad(45)

Introducing a bivector variation j = \bar{j} + \epsilon, and writing I = m l^2 \sin^2\theta_0 we have

\begin{aligned}\bar{S} + \delta S &= -\frac{1}{{2}} I \int (\bar{j}' + \epsilon')^2 dt \\ &= -\frac{1}{{2}} I \int (\bar{j}')^2 dt-\frac{1}{{2}} I \int \left( {\bar{j}}' \epsilon' + \epsilon' {\bar{j}}' \right) dt-\frac{1}{{2}} I \int (\epsilon')^2 dt\end{aligned}

The first term is just \bar{S}. Setting the variation \delta S = 0 (neglecting the quadratic \epsilon' term) and integrating by parts we have

\begin{aligned}0 &= \delta S \\ &= \int dt \left( \frac{d}{dt}\left( \frac{1}{{2}} I {\bar{j}}' \right) \epsilon + \epsilon \frac{d}{dt}\left( \frac{1}{{2}} I {\bar{j}}' \right) \right) \\ &= \int \left( \frac{d}{dt} I {\bar{j}}' \right) \dot{c} \epsilon dt \end{aligned}

With \epsilon arbitrary it appears that the solutions of the variation problem are given by

\begin{aligned}\frac{d}{dt}\left( I j' \right) = 0\end{aligned} \quad\quad\quad(46)

Has anything been lost by requiring that this is zero identically when all we had originally was the dot product of this with the variation bivector was zero? If this zero describes the solution set, then we should be able to integrate this, yielding a constant. However, following this to its logical conclusion leads to inconsistency. Integrating 46, producing a bivector constant \kappa we have

\begin{aligned}I j' = \kappa\end{aligned} \quad\quad\quad(47)

The original constraint on j was the bivector spanning the plane from the polar axis to the azimuthal unit vector in the x,y plane at angle a \phi. The \phi dependence wasn’t specified, and left encoded in the bivector representation. Without words that was

\begin{aligned}j = \mathbf{e}_3 \mathbf{e}_1 e^{\mathbf{e}_{12} \phi}\end{aligned} \quad\quad\quad(48)

Inserting back into 47 this gives

\begin{aligned}\frac{d}{dt} e^{\mathbf{e}_{12} \phi} = \frac{\mathbf{e}_1 \mathbf{e}_3 \kappa}{I}\end{aligned} \quad\quad\quad(49)

One more integration is trivially possible yielding

\begin{aligned}e^{\mathbf{e}_{12} \phi(t)} = e^{\mathbf{e}_{12} \phi_0} + \frac{\mathbf{e}_1 \mathbf{e}_3 \kappa t}{I} \end{aligned} \quad\quad\quad(50)

There are two possibilities for the grades of \kappa, one is \kappa \propto \mathbf{e}_3 \mathbf{e}_1, so that the time dependence is a scalar, and the other is \kappa \propto \mathbf{e}_3 \mathbf{e}_2 so that we have an x,y plane bivector component. Allowing for both, and separating into real and imaginary parts we have

\begin{aligned}\cos\phi - \cos\phi_0 &= \frac{1}{{I}} \kappa_r t \\ \sin\phi - \sin\phi_0 &= \frac{1}{{I}} \kappa_i t\end{aligned} \quad\quad\quad(51)

This isn’t anything close to the solution that is expected if we were to start with the scalar Lagrangian for the same problem. That Lagrangian is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} I \dot{\phi}^2\end{aligned} \quad\quad\quad(53)

and the Euler Lagrange equations give us, for a scalar constant \mu

\begin{aligned}I \dot{\phi} = \mu\end{aligned} \quad\quad\quad(54)

So the solution for \phi is just

\begin{aligned}\phi - \phi_0 = \frac{\mu t}{I}\end{aligned} \quad\quad\quad(55)

In the absence of friction this makes sense. Our angle increases monotonically, and we have circular motion with constant angular velocity. Compare this to the messier 51 derived from the bivector “solution” to the variational problem. There is definitely something wrong with the variational approach (or conclusion) when the variable is a bivector.

Does it help to include the constraints explicitly? The bivector parametrized Lagrangian with the unit bivector multiplier constraint for this system is

\begin{aligned}\mathcal{L} = -\frac{1}{{2}} I (j')^2 + \lambda (j^2 + 1)\end{aligned} \quad\quad\quad(56)

Doing the variation we get a set of two equations

\begin{aligned}- I j'' &= 2 \lambda j \\ j^2 &= -1\end{aligned} \quad\quad\quad(57)

Once again re-inserting j = \mathbf{e}_{31} e^{\mathbf{e}_{12} \phi} one gets

\begin{aligned}e^{i\phi} = A \sin(\sqrt{I/2\lambda} t) + B \cos(\sqrt{I/2\lambda} t)\end{aligned} \quad\quad\quad(59)

Setting \lambda = 1, A = i, B = 1, we have

\begin{aligned}\phi \propto t\end{aligned} \quad\quad\quad(60)

This is now consistent with the scalar Lagrangian treatment. In both cases we have the angle linearly proportional to time, and have a single parameter to adjust the angular velocity (used \mu in the scalar treatment above and have \lambda this time (although it was set to one here for convenience). The conclusion has to be that the requirement to include the multipliers for the constraints is absolutely necessary to get the right physics out of this bivector parametrized Lagrangian. The good news is that we can get the right physics out of a non-scalar treatment. What is a bit disturbing is that it was fairly difficult in the end to get from the results of the variation to a solution, and that this will not likely get any easier with a more complex system.

A second example.

The scalar expansion 82 of the kinetic term in our spherical polar Lagrangian shows a couple of other specific multivector functions we can consider the variation of.

We have considered an specific example of a Lagrangian function \mathcal{L} = f(j' \dot{c} j') without (yet) utilizing or deriving a multivector form of the Euler-Lagrange equations. Let’s consider a few more specific simple examples motivated by the expansion of the kinetic energy in 82. Lets start with

\begin{aligned}\mathcal{L} = \theta' a \dot{c} j'\end{aligned} \quad\quad\quad(61)

where a is a bivector constant, \theta a scalar, and j a bivector variable for the system. Expanding this around stationary points j = \bar{j} + \epsilon, and \theta = \bar{\theta} + \phi we have to first order

\begin{aligned}\int \delta \mathcal{L} dt &\approx \int \left( \phi' a \dot{c} \bar{j}' + \bar{\theta}' a \dot{c} \epsilon' \right) dt \\ &= \int \left( \phi \frac{d}{dt}( - a \dot{c} \bar{j}' ) + \frac{d}{dt}( - \bar{\theta}' a ) \dot{c} \epsilon \right) dt \\ \end{aligned}

In this simple system where both scalar variable \theta and bivector variable j are cyclic “coordinates”, we have solutions to the variational problem given by the pair of equations

\begin{aligned}- a \dot{c} j' &= \text{constant} \\ - \theta' a &= \text{constant}\end{aligned} \quad\quad\quad(62)

As to what, if anything, this particular Lagrangian (a fragment picked out of a real Kinetic energy function) represents physically that doesn’t matter so much, since the point of this example was to build up to treating the more general case where we are representing something physical.

A third example.

If we throw a small additional wrench into the problem above and allow a to be one of the variables our system is dependent on.

\begin{aligned}\mathcal{L}(\theta, \theta', a, a', j, j') = \theta' a \dot{c} j'\end{aligned} \quad\quad\quad(64)

It is less obvious how to do a first order Taylor expansion of this Lagrangian required for the variation around the stationary solution. If all the coordinates in the Lagrangian were scalars as in

\begin{aligned}\mathcal{L}(\theta, \theta', a, a', j, j') = \theta' a j'\end{aligned} \quad\quad\quad(65)

(dropping dot products now that these are all scalar variables), then the variation requires nothing abnormal. Suppose our stationary point has coordinates \bar{\theta}, \bar{a}, \bar{j}, with variations \alpha, \beta, \gamma that vanish at the extremes of the integral as usual.

With scalar variables our path is clear, and we just form

\begin{aligned}\delta S &= \int \delta \mathcal{L} \\ &= \int \left( \alpha \frac{\partial {\mathcal{L}}}{\partial {\theta}}+\alpha' \frac{\partial {\mathcal{L}}}{\partial {\theta'}}+\beta \frac{\partial {\mathcal{L}}}{\partial {a}}+\beta' \frac{\partial {\mathcal{L}}}{\partial {a'}}+\gamma \frac{\partial {\mathcal{L}}}{\partial {j}}+\gamma' \frac{\partial {\mathcal{L}}}{\partial {j'}}\right) \\ \end{aligned}

Here is it implied that all the partials are evaluated at the stationary points \bar{\theta}, \bar{a}, \bar{j}. Doing the integration by parts we have

\begin{aligned}\delta S&=\int \left( \alpha \left( \frac{\partial {\mathcal{L}}}{\partial {\theta}} - \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\theta'}} \right) +\beta \left( \frac{\partial {\mathcal{L}}}{\partial {a}} - \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {a'}} \right) +\gamma \left( \frac{\partial {\mathcal{L}}}{\partial {j}} - \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {j'}} \right) \right)\end{aligned}

Setting \delta S = 0 this produces the Euler-Lagrange equations for the system. For our specific Lagrangian this procedure gives

\begin{aligned}\int \delta \mathcal{L} &=\int \alpha' \bar{a} \bar{j}' + \bar{\theta}' \beta \bar{j}' + \bar{\theta}' \bar{a} \gamma' \\ &=\int \alpha \frac{d}{dt}\left( - \bar{a} \bar{j}' \right) + \bar{\theta}' \beta \bar{j}' + \frac{d}{dt}\left( -\bar{\theta}' \bar{a}\right) \gamma \\ \end{aligned}

Each time we do the first order expansion we are varying just one of the coordinates. It seems likely that the correct answer for the multivariable case will be

\begin{aligned}\int \delta \mathcal{L} &=\int \alpha \frac{d}{dt}\left( - \bar{a} \dot{c} \bar{j}' \right) + \beta \dot{c} (\bar{\theta}' \bar{j}') + \frac{d}{dt}\left( -\bar{\theta}' \bar{a}\right) \dot{c} \gamma \\ \end{aligned}

Thus the variational problem produces solutions (the coupled equations we still have to solve) of

\begin{aligned}- a \dot{c} j' &= \text{constant} \\ \theta' j' &= 0 \\ -\theta' a &= \text{constant}\end{aligned} \quad\quad\quad(66)

Multivector Euler-Lagrange equations.

Derivation.

Having considered a few specific Lagrangians dependent on multivector generalized “coordinates”, some basic comfort that it is at least possible has been gained. Let’s now move on to the general case. It is sufficient to consider Lagrangian functions that only depend on blades, since we can write any more general multivector as a sum of such blades

Write

\begin{aligned}\mathcal{L} = \mathcal{L}(X_k, X_k')\end{aligned} \quad\quad\quad(69)

where X_k \in \bigwedge^{m_k} is a blade of grade m_k. We can do a regular old chain rule expansion of this Lagrangian if coordinates are (temporarily) introduced for each of the blades X_k. For example, if X is a bivector in \mathbb{R}^{3}, we can write

\begin{aligned}X = \mathbf{e}_{12} a^{12} +\mathbf{e}_{13} a^{13} +\mathbf{e}_{23} a^{23} \end{aligned}

and then we can do a first order Taylor series expansion of \mathcal{L}(\bar{X} + \mathbf{\epsilon}) in terms of these coordinates. With \mathbf{\epsilon} = \mathbf{e}_{12}\epsilon^{12} + \mathbf{e}_{13}\epsilon^{13} + \mathbf{e}_{23}\epsilon^{23}, for this specific single variable Lagrangian we have

\begin{aligned}\delta \mathcal{L} &= \mathcal{L}(\bar{X} + \mathbf{\epsilon}) - \mathcal{L}(\bar{X}) \\ & \approx \epsilon^{12} {\left.{{\frac{\partial {\mathcal{L}}}{\partial {a^{12}}}}}\right\vert}_{{\bar{X}}}+ \epsilon^{23} {\left.{{\frac{\partial {\mathcal{L}}}{\partial {a^{23}}}}}\right\vert}_{{\bar{X}}}+ \epsilon^{13} {\left.{{\frac{\partial {\mathcal{L}}}{\partial {a^{13}}}}}\right\vert}_{{\bar{X}}}\end{aligned}

If we write the coordinate expansion of our blades X_k as

\begin{aligned}X_k = \sum_{b_k} {x_k}^{b_k} \sigma_{b_k}\end{aligned}

where b_k is some set of indexes like \{12, 23, 13\} from the \mathbb{R}^{3} bivector example, then our chain rule expansion of the Lagrangian to first order about a stationary point becomes

\begin{aligned}\delta \mathcal{L} &=\sum_{b_k} {\epsilon_k}^{b_k} {\left.{{ \frac{\partial {\mathcal{L}}}{\partial {{x_k}^{b_k}}} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}}+\sum_{b_k} ({\epsilon_k}^{b_k})' {\left.{{ \frac{\partial {\mathcal{L}}}{\partial {({x_k}^{b_k})'}} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}}\end{aligned} \quad\quad\quad(70)

Trying to write this beastie down with abstract indexing variables is a bit ugly. There’s something to be said for not trying to be general since writing this for the single variable bivector example was much clearer. However, having written down the ugly beastie, it can now be cleaned up nicely by introducing position and velocity gradient operators for each of the grades

\begin{aligned}\nabla_{X_k} &\equiv \sigma^{b_k} \frac{\partial }{\partial {{x_k}^{b_k}}} \\ \nabla_{X_k'} &\equiv \sigma^{b_k} \frac{\partial }{\partial {({x_k}^{b_k})'}}\end{aligned} \quad\quad\quad(71)

Utilizing this (assumed) orthonormal basis pair \sigma^{b_k} \dot{c} \sigma_{b_j} = {\delta^k}_j we have

\begin{aligned}\delta S &= \int \delta \mathcal{L} \\ &=\int \sum_k\mathbf{\epsilon}_k \dot{c} {\left.{{ \nabla_{X_k} \mathcal{L} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}}+{\mathbf{\epsilon}_k}' \dot{c} {\left.{{ \nabla_{X_k'} \mathcal{L} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}} \\ &=\int \sum_k\mathbf{\epsilon}_k \dot{c} \left( {\left.{{ \nabla_{X_k} \mathcal{L} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}} - \frac{d}{dt} {\left.{{ \nabla_{X_k'} \mathcal{L} }}\right\vert}_{{\bar{X}_k, \bar{X}_k'}} \right) \\ \end{aligned}

Setting \delta S = 0 we have a set of equations for each of the blade variables, and the Euler-Lagrange equations take the form

\begin{aligned}\nabla_{X_k} \mathcal{L} = \frac{d}{dt} \nabla_{X_k'} \mathcal{L} \end{aligned} \quad\quad\quad(73)

Work some examples.

For the pendulum problem we are really only interested in applying 73 for scalar and bivector variables. Let’s revisit some of the example Lagrangians, functions of bivectors, already considered and see how we can evaluate the position and velocity gradients. This generalized Euler-Lagrange formulation doesn’t do much good if we can’t work with it.

Let’s take the bivector gradient of a couple bivector functions (all the ones considered previously leading up to the multivector Euler-Lagrange equations).

\begin{aligned}f(B) &= A \dot{c} B \\ g(B) &= B^2 \end{aligned}

It is actually sufficient (in a less than 4D space) to consider only the first,

\begin{aligned}\nabla_B f &= \nabla_B \left\langle{{ A B }}\right\rangle \\ &= \sum_{a<b} \sigma^{ab} \frac{\partial }{\partial {b^{ab}}} \left\langle{{ A b^{a' b'} \sigma_{a' b'} }}\right\rangle \\ &= \sum_{a<b} \sigma^{ab} \left\langle{{ A \sigma_{ab} }}\right\rangle \\ &= \sum_{a<b} \sigma^{ab} A \dot{c} \sigma_{ab} \\ &= A\end{aligned}

For g we then have

\begin{aligned}\nabla_B g &= \nabla_B B B \\ &= \nabla_B \left\langle{{ B B}}\right\rangle \\ &= 2 B\end{aligned}

This is now enough to evaluate our bivector parameterized Lagrangian from the first example, 56, reproducing the result obtained by direct variation (as in Feynman’s lectures)

\begin{aligned}\frac{d}{dt} \nabla_{j'} \mathcal{L} &= \nabla_{j} \mathcal{L} \\ - I j'' &= \\         &= 2 \lambda j\end{aligned}

By inspection we can see that this works for the remaining two motivational examples too.

Evaluating the pendulum Euler-Lagrange equations (scalar, bivector parameterized KE).

We now have the tools to evaluate the Euler-Lagrange equations for the pendulum problem 37. Since all the d{\theta}_a/dt and dj_a/dt dependence is in the kinetic energy we can start with that. For

\begin{aligned}K = \frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\dagger Q \dot{\boldsymbol{\Theta}}\end{aligned} \quad\quad\quad(74)

We want each of \frac{\partial {K}}{\partial {\dot{\theta}_a}} and \nabla_{j_a'} K. Of these the \dot{\theta} derivative is easier, so lets start with that

\begin{aligned}\frac{\partial {K}}{\partial {\theta_a}}&=\frac{1}{{2}} \left( \frac{\partial {\dot{\boldsymbol{\Theta}}}}{\partial {\dot{\theta}_a}}^\dagger \right) Q \dot{\boldsymbol{\Theta}}+\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\dagger Q \left( \frac{\partial {\dot{\boldsymbol{\Theta}}}}{\partial {\dot{\theta}_a}} \right)  \\ \end{aligned}

Each of these are scalars and thus equal their Hermitian conjuagate. This leaves us with just one term doubled

\begin{aligned}\frac{\partial {K}}{\partial {\theta_a}}&=\dot{\boldsymbol{\Theta}}^\dagger Q \left( \frac{\partial {\dot{\boldsymbol{\Theta}}}}{\partial {\dot{\theta}_a}} \right) \\ \end{aligned}

A final evaluation of the derivatives, in block matrix form over rows r, gives us

\begin{aligned}\frac{\partial {K}}{\partial {\theta_a}}=\dot{\boldsymbol{\Theta}}^\dagger Q {\begin{bmatrix}\delta_{ar}\begin{bmatrix}1 \\ 0\end{bmatrix}\end{bmatrix}}_{r}\end{aligned} \quad\quad\quad(75)

For the bivector velocity gradients, we can do the same,

\begin{aligned}\nabla_{j_a'} K&=\frac{1}{{2}} \left( \stackrel{ \rightarrow }\nabla_{j_a'} \dot{\boldsymbol{\Theta}}^\dagger \right) Q \dot{\boldsymbol{\Theta}}+\frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\dagger Q \left( \dot{\boldsymbol{\Theta}} \stackrel{ \leftarrow }\nabla_{j_a'} \right)  \\ \end{aligned}

Only the j_a' parts of \dot{\boldsymbol{\Theta}} contribute to the velocity gradients, so we need to know how to evaluate a bivector gradient (as opposed to the square which was done earlier). Expanding in coordinates

\begin{aligned}\nabla_j j &= \frac{1}{{2}} \sum_{a b} \sigma^{ab} \frac{\partial }{\partial {j^{ab}}} \frac{1}{{2}} \sum_{a'b'} \sigma_{a' b'} j^{a' b'} \\ &= \frac{1}{{2}} \sum_{a b} \sigma^{ab} \sigma_{a b} \\ &= 3\end{aligned}

(three here because we are in \mathbb{R}^{3}, and our bivectors have three independent coordinates).

\begin{aligned}\nabla_{j_a'} K&=\frac{3}{2} {\begin{bmatrix}\delta_{ac}\begin{bmatrix}0 & -1 \end{bmatrix}\end{bmatrix}}_{c}Q \dot{\boldsymbol{\Theta}}+\frac{3}{2} \dot{\boldsymbol{\Theta}}^\dagger Q {\begin{bmatrix}\delta_{ar}\begin{bmatrix}0 \\ 1 \end{bmatrix}\end{bmatrix}}_{r}\end{aligned}

These terms are both one-by-one bivector matrixes, and negate with Hermitian conjugation, so we can again double up and eliminate one of the two terms, producing

\begin{aligned}\nabla_{j_a'} K&= 3 \dot{\boldsymbol{\Theta}}^\dagger Q {\begin{bmatrix}\delta_{ar}\begin{bmatrix}0 \\ 1 \end{bmatrix}\end{bmatrix}}_{r}\end{aligned} \quad\quad\quad(76)

Completing the Euler-Lagrange equation evaluation we have for the \theta_a coordinates

\begin{aligned}\frac{d}{dt} \left(\dot{\boldsymbol{\Theta}}^\dagger Q {\begin{bmatrix}\delta_{ar}\begin{bmatrix}1 \\ 0\end{bmatrix}\end{bmatrix}}_{r}\right)= \frac{1}{{2}} \dot{\boldsymbol{\Theta}}^\dagger \frac{\partial {Q}}{\partial {\theta_a}} \dot{\boldsymbol{\Theta}} - \frac{\partial {\phi}}{\partial {\theta_a}}\end{aligned} \quad\quad\quad(77)

plus one equation for each of the bivectors j_a

\begin{aligned} 3 \frac{d}{dt} \left(\dot{\boldsymbol{\Theta}}^\dagger Q {\begin{bmatrix}\delta_{ar}\begin{bmatrix}0 \\ 1 \end{bmatrix}\end{bmatrix}}_{r}\right)=\frac{1}{{2}} \sum_{e < f} \mathbf{e}_f \mathbf{e}_e \dot{\boldsymbol{\Theta}}^\dagger \frac{\partial {Q}}{\partial {{j_a}^{e f}}} \dot{\boldsymbol{\Theta}} + 2 \lambda_a j_a\end{aligned} \quad\quad\quad(78)

Because the bivector \mathbf{e}_f \mathbf{e}_e does not (nessessarily) commute with bivectors j_a' that are part of \dot{\boldsymbol{\Theta}} there doesn’t look like there’s much hope of assembling the left hand j_a gradients into a nice non-coordinate form. Additionally, looking back at this approach, and the troubles with producing meaningful equations of motion even in the constrained single pendulum case, it appears that a bivector parameterization of the kinetic energy is genally not a good approach. This is an unfortunate conclusion after going through all the effort to develop intuition that led to the multivector Euler-Lagrange formulation for this problem.

Oh well.

The Hermitian formulation used here should still provide a nice compact way of expressing the kinetic energy, even if we work with plain old scalar spherical polar angles \theta, and \phi. Will try this another day, since adding that to these notes will only make them that much more intractable.

Appendix calculation. A verification that the Kinetic matrix product is a real scalar.

In the kinetic term of the rather scary looking Lagrangian of 37 we have what should be a real scalar, but it is not obvious that this is the case. As a validation that nothing very bad went wrong, it seems worthwhile to do a check that this is in fact the case, expanding this out explicitly in gory detail.

One way to try this expansion is utilizing a block matrix summing over the diagonal and paired skew terms separately. That is

\begin{aligned}K &=\frac{1}{{2}}\sum_{k=1}^N\mu_k{l_k}^2{\dot{\boldsymbol{\Theta}}_k}^\dagger\begin{bmatrix}1 & - j_k e^{-j_k \theta_k} \sin\theta_k \\ j_k e^{j_k \theta_k} \sin\theta_k & \sin^2\theta_k \end{bmatrix}\dot{\boldsymbol{\Theta}}_k \\ &+\frac{1}{{2}}\sum_{a<b}\mu_bl_a l_b\left({\dot{\boldsymbol{\Theta}}_a}^\dagger\begin{bmatrix}- j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} & - j_a e^{-j_a \theta_a} \sin\theta_b \\ j_b e^{j_b \theta_b} \sin\theta_a & \sin\theta_a \sin\theta_b\end{bmatrix}\dot{\boldsymbol{\Theta}}_b +{\dot{\boldsymbol{\Theta}}_b}^\dagger\begin{bmatrix}- j_b e^{-j_b \theta_b} j_a e^{j_a\theta_a} & - j_b e^{-j_b \theta_b} \sin\theta_a \\ j_a e^{j_a \theta_a} \sin\theta_b & \sin\theta_b \sin\theta_a\end{bmatrix}\dot{\boldsymbol{\Theta}}_a\right)\end{aligned}

Examining the diagonal matrix products and expanding one of these (dropping the k suffix for tidiness), we have

\begin{aligned}{\dot{\boldsymbol{\Theta}}}^\dagger\begin{bmatrix}1 & - j e^{-j \theta} \sin\theta \\ j e^{j \theta} \sin\theta & \sin^2\theta \end{bmatrix}\dot{\boldsymbol{\Theta}} &=\dot{\theta}^2 -\sin^2\theta \left(\frac{dj}{dt}\right)^2 - \dot{\theta} \sin\theta\cos\theta \left( j \frac{dj}{dt} + \frac{dj}{dt} j \right)\end{aligned} \quad\quad\quad(79)

Since we are working in 3D this symmetric sum is twice the dot product of the bivector j with its derivative, which means that it is a scalar. We expect this to be zero though, and can observe that this is the case since j was by definition a unit bivector

\begin{aligned}j \frac{dj}{dt} + \frac{dj}{dt} j = \frac{d j^2}{dt} = \frac{d (-1)}{dt} = 0\end{aligned} \quad\quad\quad(80)

(thus j and its derivative represent orthogonal oriented planes rather like \hat{\mathbf{r}} and its derivative are orthogonal on a circle or sphere). The implication is that the diagonal subset of the kinetic energy expansion contains just

\begin{aligned}\frac{1}{{2}}\sum_{k=1}^N\mu_k{l_k}^2\left(\left(\frac{d\theta_k}{dt}\right)^2 -\sin^2\theta_k \left(\frac{d j_k}{dt}\right)^2  \right)\end{aligned} \quad\quad\quad(81)

If we are going to have any complex interaction terms then they will have to come from the off diagonal products. Expanding the first of these

\begin{aligned}{\dot{\boldsymbol{\Theta}}_a}^\dagger&\begin{bmatrix}- j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} & - j_a e^{-j_a \theta_a} \sin\theta_b \\ j_b e^{j_b \theta_b} \sin\theta_a & \sin\theta_a \sin\theta_b\end{bmatrix}\dot{\boldsymbol{\Theta}}_b \\ &=\begin{bmatrix}\theta_a' & -j_a'\end{bmatrix}\begin{bmatrix}- j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} & - j_a e^{-j_a \theta_a} \sin\theta_b \\ j_b e^{j_b \theta_b} \sin\theta_a & \sin\theta_a \sin\theta_b\end{bmatrix}\begin{bmatrix}\theta_b' \\ j_b'\end{bmatrix} \\ &=\begin{bmatrix}\theta_a' & -j_a'\end{bmatrix}\begin{bmatrix}- j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} \theta_b' - j_a e^{-j_a \theta_a} \sin\theta_b j_b' \\ j_b e^{j_b \theta_b} \sin\theta_a \theta_b' + \sin\theta_a \sin\theta_b j_b'\end{bmatrix} \\ &=- j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} \theta_a' \theta_b' - j_a e^{-j_a \theta_a} \sin\theta_b \theta_a' j_b' -j_a' j_b e^{j_b \theta_b} \sin\theta_a \theta_b' - \sin\theta_a \sin\theta_b j_a' j_b'\end{aligned}

Adding to this the a \leftrightarrow b exchanged product and rearranging yields

\begin{aligned}- \theta_a' \theta_b' (j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} + j_b e^{-j_b \theta_b} j_a e^{j_a\theta_a} ) - \sin\theta_a \sin\theta_b (j_a' j_b' + j_b' j_a') \\ - \sin\theta_b \theta_a' ( j_a e^{-j_a \theta_a} j_b' + j_b' j_a e^{j_a \theta_a} )- \sin\theta_a \theta_b' ( j_b e^{-j_b \theta_b} j_a' + j_a' j_b e^{j_b \theta_b} )\end{aligned}

Each of these multivector sums within the brackets is of the form A + \tilde{A}, a multivector plus its reverse. There can therefore be no bivector or trivector terms since they negate on reversal, and the resulting sum can have only scalar and vector grades. Of these the second term, j_a' j_b' + j_b' j_a' = 2 j_a' \dot{c} j_b' so it is unarguably a scalar as expected, but additional arguments are required to show this of the other three terms. Of these remaining three, the last two have the same form. Examining the first of these two

\begin{aligned}j_b e^{-j_b \theta_b} j_a' + j_a' j_b e^{j_b \theta_b} &=(j_b \cos\theta_b + \sin\theta_b) j_a' + j_a' (j_b \cos\theta_b - \sin\theta_b) \\ &=\cos\theta_b (j_b j_a' + j_a' j_b ) \\ &=2 \cos\theta_b (j_b \dot{c} j_a')\end{aligned}

The first term actually expands in a similarly straightforward way. The vector terms all cancel, and one is left with just

\begin{aligned}j_a e^{-j_a \theta_a} j_b e^{j_b\theta_b} + j_b e^{-j_b \theta_b} j_a e^{j_a\theta_a} = 2 \cos\theta_a\cos\theta_b j_a \dot{c} j_b - 2 \sin\theta_a \sin\theta_b \end{aligned}

Writing S_{\theta_k} = \sin\theta_k and C_{\theta_k} = \cos\theta_k (for compactness to fit things all in since the expanded result is messy), all of this KE terms can be assembled into the following explicit scalar expansion

\begin{aligned}\begin{array}{l l}K &=\frac{1}{{2}}\sum_{k=1}^N\mu_k{l_k}^2\left( (\theta_k')^2 - (S_{\theta_k} j_k')^2 \right) \\ &-\sum_{a<b}\mu_bl_a l_b\left(\theta_a'  \theta_b'  (C_{\theta_a}C_{\theta_b} (j_a \dot{c} j_b) - S_{\theta_a} S_{\theta_b} )+S_{\theta_a} S_{\theta_b} (j_a'  \dot{c} j_b' )+S_{\theta_b} C_{\theta_b} \theta_a'  (j_b \dot{c} j_a' )+S_{\theta_a} C_{\theta_a} \theta_b'  (j_a \dot{c} j_b' )\right)\end{array}\end{aligned} \quad\quad\quad(82)

Noting that all the bivector, bivector dot products are scalars really completes the desired verification. We can however, be more explicit using j_a = \mathbf{e}_{31} e^{\mathbf{e}_{12}\phi_a}, which gives after a bit of manipulation

\begin{aligned}j_a \dot{c} j_b &= - \cos(\phi_a - \phi_b) \\ (j_a')^2 &= -(\phi_a')^2 \\ j_a \dot{c} j_b' &= \phi_b' \sin(\phi_a - \phi_b) \\ j_a' \dot{c} j_b' &= -\phi_a' \phi_b' \cos(\phi_a - \phi_b)\end{aligned} \quad\quad\quad(83)

These can then be inserted back into 82 in a straightforward fashion, but it is not any more illuminating to do so.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: