# Peeter Joot's (OLD) Blog.

• ## Archives

 Adam C Scott on avoiding gdb signal noise… Ken on Scotiabank iTrade RESP …… Alan Ball on Oops. Fixing a drill hole in P… Peeter Joot's B… on Stokes theorem in Geometric… Exploring Stokes The… on Stokes theorem in Geometric…

• 320,607

## Hamiltonian notes. (a start)

Posted by peeterjoot on October 1, 2009

# Motivation

I’ve now seen Hamiltonian’s used, mostly in a Quantum context, and think that I understand at least some of the math associated with the Hamiltonian and the Hamiltonian principle. I have, however, not used either of these enough that it seems natural to do so.

Here I attempt to summarize for myself what I know about Hamiltonian’s, and work through a number of examples. Some of the examples considered will be ones already treated with the Lagrangian formalism [1].

Some notation will be invented along the way as reasonable, since I’d like to try to also relate the usual coordinate representation of the Hamiltonian, the Hamiltonian principle, and the Poisson bracket, with the bivector representation of the 2N complex configuration space introduced in [2]. (NOT YET DONE).

# Hamiltonian as a conserved quantity.

Starting with the Lagrangian formalism the Hamiltonian can be found as a conserved quantity associated with time translation when the Lagrangian has no explicit time dependence. This follows directly by considering the time derivative of the Lagrangian $\mathcal{L} = \mathcal{L}(q^i, \dot{q}^i)$.

\begin{aligned}\frac{d\mathcal{L}}{dt} &= \frac{\partial {\mathcal{L}}}{\partial {q^i}} \frac{dq^i}{dt} +\frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \frac{d\dot{q}^i}{dt} \\ &= \dot{q}^i \frac{d}{dt}\frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} +\frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \frac{d\dot{q}^i}{dt} \\ &= \frac{d}{dt} \left( \dot{q}^i \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \right) \end{aligned}

We can therefore form the difference

\begin{aligned}\frac{d}{dt} \left( \dot{q}^i \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} -\mathcal{L} \right) = 0\end{aligned} \quad\quad\quad(1)

and find that this quantity, labeled H, is a constant of motion for the system

\begin{aligned}H \equiv \dot{q}^i \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} -\mathcal{L} = \text{constant}\end{aligned} \quad\quad\quad(2)

We’ll see later that this constant is sometimes the total energy of the system.

The $\dot{q}^i$ partials of the Lagrangian are called the canonical momentum conjugate to $q^i$. Quite a mouthful, so just canonical momenta seems like a good compromise. We will write (reserving $p^i = m q^i$ for the non-canonical momenta)

\begin{aligned}P_i \equiv \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}}\end{aligned} \quad\quad\quad(3)

and note that these are the coordinates of a sort of velocity gradient of the Lagrangian. We’ve seen these canonical momenta in velocity gradient form previously where it was noted that we could write the Euler-Lagrange equations in vector form in an orthonormal reciprocal frame space as

\begin{aligned}\nabla \mathcal{L} = \frac{d}{dt} \nabla_v \mathcal{L}\end{aligned} \quad\quad\quad(4)

where $\nabla_v = e^i \partial \mathcal{L}/\partial \dot{x}^i = e^i P_i$, $\nabla = e^i \partial/\partial x^i$, and $x = e_i x^i$.

# Some syntactic sugar. In vector form.

Following Jackson [3] (section 12.1, relativistic Lorentz force Hamiltonian), this can be written in vector form if the velocity gradient, the vector sum of the momenta conjugate to the $q^i$‘s is given its own symbol $\mathbf{P}$. He writes

\begin{aligned}H = \mathbf{v} \dot{c} \mathbf{P} - \mathcal{L}\end{aligned} \quad\quad\quad(5)

This makes most sense when working in orthonormal coordinates, but can be generalized. Suppose we introduce a pair of reciprocal frame basis for the generalized position and velocity coordinates, writing as vectors in configuration space

\begin{aligned}q &= e_i q^i \\ v &= f_i \dot{q}^i \end{aligned} \quad\quad\quad(6)

Following [2] (who use this for their bivector complexification of the configuration space), we have the freedom to impose orthonormal constraints on this configuration space basis

\begin{aligned}e^i \dot{c} e_j &= {\delta^i}_j \\ f^i \dot{c} f_j &= {\delta^i}_j \\ e^i \dot{c} f_j &= {\delta^i}_j\end{aligned} \quad\quad\quad(8)

We can now define configuration space position and velocity gradients

\begin{aligned}\nabla &\equiv e^i \frac{\partial {}}{\partial {q^i}} \\ \nabla_v &\equiv f^i \frac{\partial {}}{\partial {\dot{q}^i}}\end{aligned} \quad\quad\quad(11)

so the conjugate momenta in vector form is now

\begin{aligned}P \equiv \nabla_v \mathcal{L} = f^i \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}}\end{aligned} \quad\quad\quad(13)

Our Hamiltonian takes the form

\begin{aligned}H = v \dot{c} P - \mathcal{L}\end{aligned} \quad\quad\quad(14)

# The Hamiltonian principle.

We want to take partials of 2 with respect to $P_i$ and $q^i$. In terms of the canonical momenta we want to differentiate

\begin{aligned}H \equiv \dot{q}^i P_i -\mathcal{L}(q^i, \dot{q}^i, t)\end{aligned} \quad\quad\quad(15)

for the $P_i$ partial we have

\begin{aligned}\frac{\partial {H}}{\partial {P_i}} = \dot{q}^i\end{aligned}

and for the $q^i$ partial

\begin{aligned}\frac{\partial {H}}{\partial {q^i}} &= -\frac{\partial {\mathcal{L}}}{\partial {q^i}} \\ &= - \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \end{aligned}

These two results taken together form what I believe is called the Hamiltonian principle

\begin{aligned}\frac{\partial {H}}{\partial {P_i}} &= \dot{q}^i \\ \frac{\partial {H}}{\partial {q^i}} &= - \dot{P}_i \\ P_i &= \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \end{aligned} \quad\quad\quad(16)

A set of 2N first order equations equivalent to the second order Euler-Lagrange equations. These appear to follow straight from the definitions. Given that I’m curious why the more complex method of derivation is chosen in [4]. There the total differential of the Hamiltonian is computed

\begin{aligned}dH &= \dot{q}^i dP_i + d\dot{q}^i P_i - dq^i \frac{\partial {\mathcal{L}}}{\partial {q^i}}- d \dot{q}^i \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}}- dt \frac{\partial {\mathcal{L}}}{\partial {t}} \\ &= \dot{q}^i dP_i + d\dot{q}^i \left( P_i - \frac{\partial {\mathcal{L}}}{\partial {\dot{q}^i}} \right)- dq^i \frac{\partial {\mathcal{L}}}{\partial {q^i}}- dt \frac{\partial {\mathcal{L}}}{\partial {t}} \\ &= \dot{q}^i dP_i - dq^i \underbrace{\frac{\partial {\mathcal{L}}}{\partial {q^i}}}_{= d/dt P_i}- dt \frac{\partial {\mathcal{L}}}{\partial {t}} \\ \end{aligned}

A term by term comparison to the total differential written out explicitly

\begin{aligned}dH &= \frac{\partial {H}}{\partial {q^i}} d q^i+\frac{\partial {H}}{\partial {P_i}} d P_i+\frac{\partial {H}}{\partial {t}} dt\end{aligned} \quad\quad\quad(19)

allows the Hamiltonian equations to be picked off.

\begin{aligned}\frac{\partial {H}}{\partial {P_i}} &= \dot{q}^i \\ \frac{\partial {H}}{\partial {q^i}} &= - \dot{P}_i \\ \frac{\partial {H}}{\partial {t}} &= - \frac{\partial {\mathcal{L}}}{\partial {t}} \end{aligned} \quad\quad\quad(20)

I guess that isn’t that much more complicated and it does yield a relation between the Hamiltonian and Lagrangian time derivatives.

# Examples.

Now, that’s just about the most abstract way we can start things off isn’t it? Getting some initial feel for this constant of motion can be had by considering a sequence of Lagrangians, starting with the very simplest.

## Force free motion

Our very simplest Lagrangian is that of one dimensional purely kinetic motion

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m v^2 = \frac{1}{{2}} m \dot{x}^2\end{aligned} \quad\quad\quad(23)

Our Hamiltonian is in this case just

\begin{aligned}H = \dot{x} m \dot{x} - \frac{1}{{2}} m \dot{x} = \frac{1}{{2}} m v^2\end{aligned} \quad\quad\quad(24)

The Hamiltonian is just the kinetic energy. The canonical momentum in this case is also equal to the momentum, so eliminating $v$ to apply the Hamiltonian equations we have

\begin{aligned}H = \frac{1}{{2m}} p^2\end{aligned} \quad\quad\quad(25)

We have then

\begin{aligned}\frac{\partial {H}}{\partial {p}} &= \frac{p}{m} = \dot{x} \\ \frac{\partial {H}}{\partial {x}} &= 0 = -\dot{p} \end{aligned}

Just for fun we can put this simple linear system in matrix form

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p \\ x\end{bmatrix}=\frac{1}{m}\begin{bmatrix}0 & 0 \\ 1 & 0\end{bmatrix}\begin{bmatrix}p \\ x\end{bmatrix}\end{aligned} \quad\quad\quad(26)

A linear system of this form $y' = A y$ can be solved by exponentiation with solution

\begin{aligned}y = e^{A t} y_0\end{aligned} \quad\quad\quad(27)

In this case our matrix is nilpotent degree 2 so we can exponentiate only requiring up to the first order power

\begin{aligned}e^{A t} = I + A t\end{aligned} \quad\quad\quad(28)

specifically

\begin{aligned}\begin{bmatrix}p \\ x\end{bmatrix}=\begin{bmatrix}1 & 0 \\ \frac{t}{m} & 1\end{bmatrix}\begin{bmatrix}p_0 \\ x_0\end{bmatrix}\end{aligned} \quad\quad\quad(29)

Written out in full this is just

\begin{aligned}p &= p_0 \\ x &= \frac{p_0}{m} t + x_0\end{aligned} \quad\quad\quad(30)

Since the canonical momentum is the regular momentum $p = m v$ in this case, we have the usual constant rate change of position $x = v_0 t + x_0$ that we could have gotten in many easier ways. I’d hazard a guess that any single variable Lagrangian that is at most quadratic in position or velocity will yield a linear system.

The generalization of this Hamiltonian to three dimensions is straightforward, and we get

\begin{aligned}H &= \frac{1}{m} \mathbf{p}^2 \end{aligned} \quad\quad\quad(32)

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p_x \\ x \\ p_y \\ y \\ p_z \\ z \\ \end{bmatrix}=\frac{1}{m}\begin{bmatrix}0 & 0 & & & & \\ 1 & 0 & & & & \\ & & 0 & 0 & & \\ & & 1 & 0 & & \\ & & & & 0 & 0 \\ & & & & 1 & 0 \\ \end{bmatrix}\begin{bmatrix}p_x \\ x \\ p_y \\ y \\ p_z \\ z \\ \end{bmatrix}\end{aligned} \quad\quad\quad(33)

Since there is no coupling (nilpotent matrices down the diagonal) between the coordinates this can be treated as three independent sets of equations of the form 26, and we have

\begin{aligned}p_i(t) &= p_i(0) \\ x_i(t) &= \frac{p_i(0)}{m} t + x_i(0)\end{aligned} \quad\quad\quad(34)

Or just

\begin{aligned}\mathbf{p}(t) &= \mathbf{p}(0) \\ \mathbf{x}(t) &= \frac{\mathbf{p}(0)}{m} t + \mathbf{x}(0)\end{aligned} \quad\quad\quad(36)

## Linear potential (surface gravitation)

For the gravitational force $F = - m g \hat{\mathbf{z}} = - \boldsymbol{\nabla} \phi$, we have $\phi = m g z$, and a Lagrangian of

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m \mathbf{v}^2 - \phi = \frac{1}{{2}} m \mathbf{v}^2 - m g z\end{aligned} \quad\quad\quad(38)

Without velocity dependence the canonical momentum is the momentum $m \mathbf{v}$, and our Hamiltonian is

\begin{aligned}H = \frac{1}{{2 m}} \mathbf{p}^2 + m g z\end{aligned} \quad\quad\quad(39)

The Hamiltonian equations are

\begin{aligned}\frac{\partial {H}}{\partial {p_i}} &= \dot{x}_i = \frac{1}{m} p_i \\ \sigma_i \frac{\partial {H}}{\partial {x_i}} &= -\sigma_i \dot{p}_i = \begin{bmatrix}0 \\ 0 \\ m g \end{bmatrix}\end{aligned} \quad\quad\quad(40)

In matrix form we have

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p_x \\ x \\ p_y \\ y \\ p_z \\ z \\ \end{bmatrix}=\frac{1}{m}\begin{bmatrix}0 & 0 & & & & \\ 1 & 0 & & & & \\ & & 0 & 0 & & \\ & & 1 & 0 & & \\ & & & & 0 & 0 \\ & & & & 1 & 0 \\ \end{bmatrix}\begin{bmatrix}p_x \\ x \\ p_y \\ y \\ p_z \\ z \\ \end{bmatrix}+\begin{bmatrix}0 \\ 0 \\ 0 \\ 0 \\ -m g \\ 0 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(42)

So our problem is now reduced to solving a linear system of the form

\begin{aligned}y' = A y + b\end{aligned} \quad\quad\quad(43)

That extra little term $b$ throws a wrench into things and I’m no longer sure how to integrate by inspection. What can be noted is that we really only have to consider the $z$ components since we’ve solved the problem for the $x$ and $y$ coordinates in the force free case. That leaves

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p_z \\ z \\ \end{bmatrix}=\frac{1}{m}\begin{bmatrix}0 & 0 \\ 1 & 0 \\ \end{bmatrix}\begin{bmatrix}p_z \\ z \\ \end{bmatrix}+\begin{bmatrix}-m g \\ 0 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(44)

Is there any reason that we have to solve in matrix form? Except for a coolness factor, not really, and we can integrate each equation directly. For the momentum equation we have

\begin{aligned}p_z = - m g t + p_z(0)\end{aligned} \quad\quad\quad(45)

This can be substituted into the position equation for

\begin{aligned}\dot{z} = \frac{1}{m} (p_z(0) - m g t)\end{aligned} \quad\quad\quad(46)

Direct integration is now possible for the final solution

\begin{aligned}z &= \frac{1}{m} (p_z(0) t - m g t^2/2) + z_0 \\ &= \frac{p_z(0)}{m} t - \frac{g}{2} t^2 + z_0 \end{aligned}

Again something that we could have gotten in many easier ways. Using the result we see that the solution to 44 in matrix form, again with $A = \frac{1}{m}\begin{bmatrix}0 & 0 \\ 1 & 0\end{bmatrix}$ is

\begin{aligned}\begin{bmatrix}p_z \\ z \\ \end{bmatrix}= e^{At} \begin{bmatrix}p_z(0) \\ z(0) \\ \end{bmatrix}- m g \begin{bmatrix}t \\ \frac{1}{{2m}} t^2\end{bmatrix}\end{aligned} \quad\quad\quad(47)

I thought if I wrote this out how to solve 43 may be more obvious, but that path is still unclear. If $A$ were invertible, which it isn’t, then writing $b = A c$ would allow for a change of variables. Does this matter for consideration of a physical problem. Not really, so I’ll fight the urge to play with the math for a while and perhaps revisit this later separately.

## Harmonic oscillator (spring potential)

Like the free particle, the harmonic oscillator is very tractable in a phase space representation. For a restoring force $F = - k x \hat{\mathbf{x}} = -\boldsymbol{\nabla} \phi$, we have $\phi = k x^2/2$, and a Lagrangian of

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m \mathbf{v}^2 - \frac{1}{{2}} k \mathbf{x}^2\end{aligned} \quad\quad\quad(48)

Our Hamiltonian is again just the total energy

\begin{aligned}H = \frac{1}{{2m}} \mathbf{p}^2 + \frac{1}{{2}} k \mathbf{x}^2\end{aligned} \quad\quad\quad(49)

Evaluating the Hamiltonian equations we have

\begin{aligned}\frac{\partial {H}}{\partial {p_i}} &= \dot{x_i} = p_i/m \\ \frac{\partial {H}}{\partial {x_i}} &= -\dot{p_i} = k x_i\end{aligned} \quad\quad\quad(50)

Considering just the $x$ dimension (the others have the free particle behavior), our matrix phase space representation is

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p \\ x \\ \end{bmatrix}=\begin{bmatrix}0 & - k \\ 1/m & 0 \\ \end{bmatrix}\begin{bmatrix}p \\ x \\ \end{bmatrix}\end{aligned} \quad\quad\quad(52)

So with

\begin{aligned}A = \begin{bmatrix}0 & - k \\ 1/m & 0 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(53)

Our solution is

\begin{aligned}\begin{bmatrix}p \\ x \\ \end{bmatrix}=e^{At}\begin{bmatrix}p_0 \\ x_0 \\ \end{bmatrix}\end{aligned} \quad\quad\quad(54)

The stateful nature of the phase space solution is interesting. Provided we can make a simultaneous measurement of position and momentum, this initial state determines a next position and momentum state at a new time $t = t_0 + \Delta t_1$, and we have a trajectory through phase space of discrete transitions from one state to another

\begin{aligned}{\begin{bmatrix}p \\ x \\ \end{bmatrix}}_{i+1}=e^{A \Delta t_{i+1}}{\begin{bmatrix}p \\ x \\ \end{bmatrix}}_i\end{aligned} \quad\quad\quad(55)

Or

\begin{aligned}{\begin{bmatrix}p \\ x \\ \end{bmatrix}}_{i+1}=e^{A \Delta t_{i+1}} e^{A \Delta t_{i}} \dot{c}s e^{A \Delta t_1}{\begin{bmatrix}p \\ x \\ \end{bmatrix}}_0\end{aligned} \quad\quad\quad(56)

As for solving the system, we require again the exponential of our matrix. This matrix being antisymmetric, has complex eigenvalues and again cannot be exponentiated easily by diagonalization. However, this antisymmetric matrix is very much like the complex imaginary and its square is a negative scalar multiple of identity, so we can proceed directly forming the power series

\begin{aligned}A^2 =\begin{bmatrix}0 & - k \\ 1/m & 0 \\ \end{bmatrix}\begin{bmatrix}0 & - k \\ 1/m & 0 \\ \end{bmatrix}=-\frac{k}{m} I\end{aligned} \quad\quad\quad(57)

The first few powers are

\begin{aligned}A^2 &= -\frac{k}{m} I \\ A^3 &= -\frac{k}{m} A \\ A^4 &= \left(\frac{k}{m}\right)^2 I \\ A^5 &= \left(\frac{k}{m}\right)^2 A\end{aligned} \quad\quad\quad(58)

So exponentiating we can collect cosine and sine terms

\begin{aligned}e^{At} &= I \left( 1 - \frac{k}{m} \frac{t^2}{2!} + \left( \frac{k}{m} \right)^2 \frac{t^4}{4!} + \dot{c}s \right) + A\sqrt{\frac{m}{k}} \left( \sqrt{\frac{k}{m}} - \left(\sqrt{\frac{k}{m}}\right)^3 \frac{t^3}{3!} + \left(\sqrt{\frac{k}{m}}\right)^5 \frac{t^5}{5!} \right) \\ &= I \cos\left(\sqrt{\frac{k}{m}} t\right) + \sqrt{\frac{m}{k}} A \sin\left(\sqrt{\frac{k}{m}} t\right)\end{aligned}

As a check it is readily verified that this satisfies the desired $d(e^{At})/dt = A e^{At}$ property.

The full solution in phase space representation is therefore

\begin{aligned}\begin{bmatrix}p \\ x \\ \end{bmatrix}=\begin{bmatrix}p_0 \\ x_0 \\ \end{bmatrix}\cos\left(\sqrt{\frac{k}{m}} t\right) + \sqrt{\frac{m}{k}}\begin{bmatrix}-k x_0 \\ p_0/m \\ \end{bmatrix}\sin\left(\sqrt{\frac{k}{m}} t\right)\end{aligned} \quad\quad\quad(62)

Written out separately this is clearer

\begin{aligned}p &= p_0 \cos\left(\sqrt{\frac{k}{m}} t\right) - \sqrt{\frac{m}{k}} k x_0 \sin\left(\sqrt{\frac{k}{m}} t\right) \\ x &= x_0 \cos\left(\sqrt{\frac{k}{m}} t\right) + \sqrt{\frac{m}{k}} \frac{p_0}{m} \sin\left(\sqrt{\frac{k}{m}} t\right)\end{aligned} \quad\quad\quad(63)

One can readily verify that $m \dot{x} = p$, and $m \dot{d}{x} = -k x$ as expected.

Let’s pause before leaving the harmonic oscillator to see if 63 seems to make sense. Consider the position solution. With only initial position and no initial velocity $p_0/m$ we have oscillation that has no dependence on the mass or spring constant. This is the unmoving mass about to be let go at the end of a spring case, and since we have no damping force the magnitude of the oscillation is exactly the initial position of the mass. If the instantaneous velocity is measured at position zero, it makes sense in this case that the oscillation amplitude does depend on both the mass and the spring constant. The stronger the spring ($k$), the bigger the oscillation, and the smaller the mass, the bigger the oscillation.

It is definitely no easier to work with the phase space formulation than just solving the second order system directly. The fact that we have a linear system to solve, at least in this particular case is kind of nice. Perhaps this methodology can be helpful considering linear approximation solutions in a neighborhood of some phase space point for more complex non-linear systems.

## Force free system dependent on only differences.

In gravitational and electrostatic problems are forces are all functions of only the difference in positions of the particles. Lets look at how the purely kinetic Lagrangian and Hamiltonian change when one or more of the vector positions is reexpressed in terms of a difference in position change of variables. In the force free case this is primarily a task of rewriting the Hamiltonian in terms of the conjugate momenta after such a change of variables.

The very simplest case is the two particle single dimensional Kinetic Lagrangian,

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_1 {\dot{r}_1}^2 + \frac{1}{{2}} m_2 {\dot{r}_2}^2 \end{aligned} \quad\quad\quad(65)

With a change of variables

\begin{aligned}x &= r_2 - r_1 \\ y &= r_2\end{aligned} \quad\quad\quad(66)

and elimination of $r_1$, and $r_2$ we have

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_1 (\dot{y} - \dot{x})^2 + \frac{1}{{2}} m_2 {\dot{y}}^2 \end{aligned} \quad\quad\quad(68)

We now need the conjugate momenta in terms of $\dot{x}$ and $\dot{y}$. These are

\begin{aligned}P_x &= \frac{\partial {\mathcal{L}}}{\partial {\dot{x}}} = -m_1 (\dot{y} - \dot{x}) \\ P_y &= \frac{\partial {\mathcal{L}}}{\partial {\dot{y}}} = m_1 (\dot{y} - \dot{x}) + m_2 \dot{y}\end{aligned} \quad\quad\quad(69)

We must now rewrite the Lagrangian in terms of $P_x$ and $P_y$, essentially requiring the inversion of this which amounts to the inversion of the two by two linear system of 69. That is

\begin{aligned}\begin{bmatrix}\dot{x} \\ \dot{y}\end{bmatrix}={\begin{bmatrix}m_1 & -m_1 \\ -m_1 & (m_1 + m_2) \\ \end{bmatrix}}^{-1}\begin{bmatrix}P_x \\ P_y \end{bmatrix}\end{aligned} \quad\quad\quad(71)

This is

\begin{aligned}\begin{bmatrix}\dot{x} \\ \dot{y}\end{bmatrix}=\frac{1}{{m_1}}\frac{1}{{m_2}}\begin{bmatrix}m_1 + m_2 & m_1 \\ m_1 & m_1 \\ \end{bmatrix}\begin{bmatrix}P_x \\ P_y \end{bmatrix}\end{aligned} \quad\quad\quad(72)

Of these only $\dot{y}$ and $\dot{y} - \dot{x}$ are of interest and after a bit of manipulation we find

\begin{aligned}\dot{y} &= \frac{1}{{m_2}}(P_x + P_y) \\ \dot{x} &= \frac{1}{{m_1}}\frac{1}{{m_2}}((m_1 + m_2)P_x + m_1 P_y)\end{aligned} \quad\quad\quad(73)

From this we find the Lagrangian in terms of the conjugate momenta

\begin{aligned}\mathcal{L} = \frac{1}{{2 m_1}} {P_x}^2 + \frac{1}{{2 m_2}} (P_x + P_y)^2\end{aligned} \quad\quad\quad(75)

A quick check shows that $P_x + P_y = m_2 \dot{r}_2$, and $P_x = -m_1 \dot{r}_1$, so we have agreement with the original Lagrangian. Generalizing to the three dimensional case is straightforward, and we have

\begin{aligned}\mathcal{L} &= \frac{1}{{2}} m_1 {\mathbf{\dot{r}}_1}^2 + \frac{1}{{2}} m_2 {\mathbf{\dot{r}}_2}^2 - \phi(\mathbf{x}_1 - \mathbf{x}_2) \end{aligned} \quad\quad\quad(76)

With

\begin{aligned}\mathbf{x} &= \mathbf{x}_1 - \mathbf{x}_2 \\ \mathbf{y} &= \mathbf{x}_2\end{aligned} \quad\quad\quad(77)

The 3D generalization of the above (following by adding indexes then summing) becomes

\begin{aligned}\mathbf{P}_x &= \sigma_j \frac{\partial {\mathcal{L}}}{\partial {\dot{x}^j}} = -m_1 (\mathbf{\dot{y}} - \mathbf{\dot{x}}) \\ \mathbf{P}_y &= \sigma_j \frac{\partial {\mathcal{L}}}{\partial {\dot{y}^j}} = m_1 (\mathbf{\dot{y}} - \mathbf{\dot{x}}) + m_2 \mathbf{\dot{y}}\end{aligned} \quad\quad\quad(79)

\begin{aligned}\mathcal{L} &= \frac{1}{{2 m_1}} {\mathbf{P}_x}^2 + \frac{1}{{2 m_2}} (\mathbf{P}_x + \mathbf{P}_y)^2 - \phi(\mathbf{x}) \\ H &= \frac{1}{{2 m_1}} {\mathbf{P}_x}^2 + \frac{1}{{2 m_2}} (\mathbf{P}_x + \mathbf{P}_y)^2 + \phi(\mathbf{x})\end{aligned} \quad\quad\quad(81)

Finally, evaluation of the Hamiltonian equations we have

\begin{aligned}\sigma_j \frac{\partial {H}}{\partial {P_x^j}} &= \mathbf{\dot{x}} \\ &= \sigma_j \left( \frac{1}{{m_1}} P_x^j + \frac{1}{{m_2}} (P_x^j + P_y^j) \right) \\ &= \frac{1}{{m_1}} \mathbf{P}_x + \frac{1}{{m_2}} (\mathbf{P}_x + \mathbf{P}_y) \end{aligned}

\begin{aligned}\sigma_j \frac{\partial {H}}{\partial {P_y^j}} &= \mathbf{\dot{y}} \\ &= \sigma_j \frac{1}{{m_2}} (P_x^j + P_y^j) \\ &= \frac{1}{{m_2}} (\mathbf{P}_x + \mathbf{P}_y) \end{aligned}

\begin{aligned}\sigma_j \frac{\partial {H}}{\partial {x^j}} &= -\mathbf{\dot{P}}_x \\ &= -\sigma_j \frac{\partial {\phi}}{\partial {x^j}} \\ &= -\boldsymbol{\nabla}_{\mathbf{x}} \phi(\mathbf{x})\end{aligned}

\begin{aligned}\sigma_j \frac{\partial {H}}{\partial {y^j}} &= -\mathbf{\dot{P}}_y \\ &= -\sigma_j \frac{\partial {\phi}}{\partial {y^j}} \\ &= 0\end{aligned}

Summarizing we have four first order equations

\begin{aligned}\mathbf{\dot{x}} &= \left( \frac{1}{{m_1}} + \frac{1}{{m_2}} \right) \mathbf{P}_x + \frac{1}{{m_2}} \mathbf{P}_y \\ \mathbf{\dot{y}} &= \frac{1}{{m_2}} (\mathbf{P}_x + \mathbf{P}_y) \\ \mathbf{\dot{P}}_x &= \boldsymbol{\nabla}_{\mathbf{x}} \phi(\mathbf{x}) \\ \mathbf{\dot{P}}_y &= 0\end{aligned} \quad\quad\quad(83)

FIXME: what would we get if using the center of mass position as one of the variables. A parametrization with three vector variables should also still work, even if it includes additional redundancy.

## Gravitational potential

Next I’d like to consider a two particle gravitational interaction. However, to start we need the Lagrangian, and what should the potential term be in a two particle gravitational Lagrangian? I’d guess something with a $1/x$ form, but do we need one potential term for each mass or something interrelated? Whatever the Lagrangian is, we want it to produce the pair of force relationships

\begin{aligned}\text{Force on 2} &= - G m_1 m_2 \frac{(\mathbf{r}_2 - \mathbf{r}_1)}{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}} \\ \text{Force on 1} &= G m_1 m_2 \frac{(\mathbf{r}_2 - \mathbf{r}_1)}{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}}\end{aligned} \quad\quad\quad(87)

Guessing that the Lagrangian has a single term for the interaction potential

\begin{aligned}\phi_{21} &= \kappa \frac{1}{{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}}} \end{aligned} \quad\quad\quad(89)

so that we have

\begin{aligned}\mathcal{L} &= \frac{1}{{2}} m {\mathbf{v}_1}^2 +\frac{1}{{2}} m {\mathbf{v}_2}^2 - \phi_{21}\end{aligned} \quad\quad\quad(90)

We can evaluate the Euler-Lagrange equations and see if the result is consistent with the Newtonian force laws of 87. Suppose we write the coordinates of $\mathbf{r}_i$ as ${x^k}_i$. There are then six Euler-Lagrange equations

\begin{aligned}\frac{\partial {\mathcal{L}}}{\partial {{x^j}_i}} &= \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {{\dot{x}^j}_i}} \\ -\frac{\partial {\phi_{21}}}{\partial {{x^j}_i}} &= m_i {\dot{d}{x}^j}_i\end{aligned}

Evaluating the potential derivatives separately. Consider the $i=2$ derivative

\begin{aligned}\frac{\partial {\phi_{21}}}{\partial {{x^j}_2}}&= \kappa \frac{\partial {}}{\partial {{x^j}_2}} \left( \sum_k ({x^k}_2 - {x^k}_1)^2 \right)^{-1/2} \\ &=-\kappa \frac{1}{{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}^3}} \sum_k ({x^k}_2 - {x^k}_1) \frac{\partial {}}{\partial {{x^j}_2}} ({x^k}_2 - {x^k}_1) \\ &=-\kappa \frac{1}{{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}^3}} ({x^j}_2 - {x^j}_1)\end{aligned}

Therefore the final result of the Euler-Lagrange equations is

\begin{aligned}\kappa \frac{1}{{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}^3}} ({x^j}_2 - {x^j}_1) &= m_2 {\dot{d}{x}^j}_2 \\ -\kappa \frac{1}{{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}^3}} ({x^j}_2 - {x^j}_1) &= m_1 {\dot{d}{x}^j}_1\end{aligned} \quad\quad\quad(91)

which confirms the Lagrangian and potential guess and fixes the constant $\kappa = - G m_1 m_2$. With the sign fixed, our potential, Lagrangian, and Hamiltonian are respectively

\begin{aligned}\phi_{21} &= -\frac{G m_2 m_1 }{{\left\lvert{\mathbf{r}_2 - \mathbf{r}_1}\right\rvert}} \\ \mathcal{L} &= \frac{1}{{2}} m_1 {\mathbf{v}_1}^2 +\frac{1}{{2}} m_2 {\mathbf{v}_2}^2 - \phi_{21} \\ H &= \frac{1}{{2 m_1}} {\mathbf{p}_1}^2 +\frac{1}{{2 m_2}} {\mathbf{p}_2}^2 + \phi_{21}\end{aligned} \quad\quad\quad(93)

There is however an undesirable asymmetry to this expression, in particular a formulation that extends to multiple particles seems desirable. Let’s write instead a slight variation

\begin{aligned}\phi_{ij} = - \frac{G m_i m_j}{{\left\lvert{\mathbf{r}_i - \mathbf{r}_j}\right\rvert}}\end{aligned} \quad\quad\quad(96)

and form a scaled by two double summation over all pairs of potentials

\begin{aligned}\mathcal{L} &= \sum_i \frac{1}{{2}} m_i {\mathbf{v}_i}^2 - \frac{1}{{2}} \sum_{i \ne j} \phi_{ij}\end{aligned} \quad\quad\quad(97)

Having established what seems like an appropriate form for the Lagrangian, we can write the Hamiltonian for the multiparticle gravitational interaction by inspection

\begin{aligned}H = \sum_i \frac{1}{{2 m_i}} {\mathbf{p}_i}^2 + \frac{1}{{2}} \sum_{i \ne j} \phi_{ij}\end{aligned} \quad\quad\quad(98)

This leaves us finally in position to evaluate the Hamiltonian equations, but the result of doing so is rudely nothing more than the Newtonian equations in coordinate form. We get, for the $k$th component of the $i$th particle

\begin{aligned}\frac{\partial {H}}{\partial {{p^k}_i}} = {\dot{x}^k}_i = \frac{1}{{m_i}} {p^k}_i \end{aligned} \quad\quad\quad(99)

\begin{aligned}\frac{\partial {H}}{\partial {{x^k}_i}} = -{\dot{p}^k}_i = G \sum_{j \ne i} m_i m_j \frac{{x^k}_i -{x^k}_j}{{\left\lvert{\mathbf{r}_i - \mathbf{r}_j}\right\rvert}^3}\end{aligned} \quad\quad\quad(100)

The state space vector for this system of equations is brutally ugly, and could be put into the following form for example

\begin{aligned}\mathbf{z} =\begin{bmatrix}{p^1}_1 \\ {p^2}_1 \\ {p^3}_1 \\ {x^1}_1 \\ {x^2}_1 \\ {x^3}_1 \\ {p^1}_2 \\ {p^2}_2 \\ {p^3}_2 \\ {x^1}_2 \\ \dot{v}s \\ \end{bmatrix}\end{aligned} \quad\quad\quad(101)

Where the Hamiltonian equations take the form of a non-linear function on such state space vectors
We have a somewhat sparse equation of the form

\begin{aligned}\frac{d\mathbf{z}}{dt} = A(\mathbf{z})\end{aligned} \quad\quad\quad(102)

One thing that is possible in such a representation is calculating the first order approximate change in position and momentum moving from one time to a small time later

\begin{aligned}\mathbf{z}(t_0 + \Delta t) = \mathbf{z}(t_0) + A(\mathbf{z}(t_0)) \Delta t\end{aligned} \quad\quad\quad(103)

One could conceivably calculate the trajectories in phase space using such increments, and if a small enough time increment is used this can be thought of as solving the gravitational system. I recall that Feynman did something like this in his lectures, but set up the problem in a more computationally efficient form (it definitely didn’t have the redundancy built into the Hamiltonian equations).

FIXME: should be able to solve this for an arbitrary $\Delta t$ later time if this was extended to the higher order terms. Need something like the $e^{z \dot{c} \nabla}$ chain rule expansion. Think this through. Will be a little different since we are already starting with the first order contribution.

What does this system of equations look like with a reduction of order through center of mass change of variables?

TODO.

TODO.

TODO.

TODO.

TODO.

TODO.

TODO.

TODO.

# References

[1] Peeter Joot. Solutions to David Tong’s mf1 Lagrangian problems [online]. http://sites.google.com/site/peeterjoot/geometric-algebra/tong_mf1.pdf.

[2] C. Doran and A.N. Lasenby. Geometric algebra for physicists. Cambridge University Press New York, Cambridge, UK, 1st edition, 2003.

[3] JD Jackson. Classical Electrodynamics Wiley. 2nd edition, 1975.

[4] H. Goldstein. Classical mechanics. Cambridge: Addison-Wesley Press, Inc, 1st edition, 1951.