Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Hamiltonian treatment of the rigid pendulum problem.

Posted by peeterjoot on October 7, 2009

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

The bob speed for a stiff rod of length l is (l \dot{\theta})^2, and our potential is m g h = m g l (1 -\cos\theta). The Lagrangian is therefore

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m l^2 \dot{\theta}^2 - mg l (1 -\cos\theta) \\ \end{aligned} \quad\quad\quad(121)

The constant m g l term can be dropped, and our canonical momentum conjugate to \dot{\theta} is p_\theta = m l^2 \dot{\theta}, so our Hamiltonian is

\begin{aligned}H = \frac{1}{{2 m l^2}} {p_\theta}^2 - mg l \cos\theta \\ \end{aligned} \quad\quad\quad(123)

We can now compute the Hamiltonian equations

\begin{aligned}\frac{\partial {H}}{\partial {p_\theta}}  &= \dot{\theta}         = \frac{1}{{ m l^2 }} p_\theta \\ \frac{\partial {H}}{\partial {q}}         &= -\dot{p}_\theta   = m g l \sin\theta\end{aligned} \quad\quad\quad(125)

Only in the neighborhood of a particular angle can we write this in matrix form. Suppose we expand this around \theta = \theta_0 + \alpha. The sine is then

\begin{aligned}\sin\theta \approx \sin\theta_0 + \cos\theta_0 \alpha\end{aligned} \quad\quad\quad(127)

The linear approximation of the Hamiltonian equations after a change of variables become

\begin{aligned}\frac{d}{dt}\begin{bmatrix}p_\theta \\ \alpha\end{bmatrix}=\begin{bmatrix}0 & -m g l \cos\theta_0 \\ 1/ m l^2 & 0\end{bmatrix}\begin{bmatrix}p_\theta \\ \alpha\end{bmatrix}-m g l \sin\theta_0\begin{bmatrix}1 \\ 0\end{bmatrix}\end{aligned} \quad\quad\quad(128)

A change of variables that scales the factors in the matrix to have equal magnitude and equivalent dimensions is helpful. Writing

\begin{aligned}\begin{bmatrix}p_\theta \\ \alpha\end{bmatrix}=\begin{bmatrix}a & 0 \\ 0 & 1\end{bmatrix}\mathbf{z}\end{aligned} \quad\quad\quad(129)

one finds

\begin{aligned}\frac{d\mathbf{z}}{dt}&=\begin{bmatrix}0 & -m g l \cos\theta_0/a \\ a/ m l^2 & 0\end{bmatrix}\mathbf{z}-\frac{m g l \sin\theta_0 }{a}\begin{bmatrix}1 \\ 0\end{bmatrix} \\ \end{aligned} \quad\quad\quad(130)

To tidy this up, we want

\begin{aligned}{\left\lvert{\frac{a}{m l^2}}\right\rvert} = {\left\lvert{\frac{m g l \cos\theta_0}{a}}\right\rvert}\end{aligned} \quad\quad\quad(132)

Or

\begin{aligned}a = m l^2 \sqrt{\frac{g}{l} {\left\lvert{\cos\theta_0}\right\rvert}}\end{aligned} \quad\quad\quad(133)

The result of applying this scaling is quite different above and below the horizontal due to the sign difference in the cosine. Below the horizontal where \theta_0 \in (-\pi/2, \pi/2) we get

\begin{aligned}\frac{d\mathbf{z}}{dt}&=\sqrt{\frac{g}{l} \cos\theta_0}\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\mathbf{z}-\sin\theta_0 \sqrt{\frac{g}{l \cos\theta_0}}\begin{bmatrix}1 \\ 0\end{bmatrix}\end{aligned} \quad\quad\quad(134)

and above the horizontal where \theta_0 \in (\pi/2, 3\pi/2) we get

\begin{aligned}\frac{d\mathbf{z}}{dt}&=\sqrt{-\frac{g}{l} \cos\theta_0}\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}\mathbf{z}-\sin\theta_0 \sqrt{-\frac{g}{l \cos\theta_0}}\begin{bmatrix}1 \\ 0\end{bmatrix}\end{aligned} \quad\quad\quad(135)

Since (\begin{smallmatrix} 0 & -1 \\  1 & 0 \end{smallmatrix}) has the characteristics of an imaginary number (squaring to the negative of the identity) the homogeneous part of the solution for the change of the phase space vector in the vicinity of any initial angle in the lower half plane is trigonometric. Similarly the solutions are necessarily hyperbolic in the upper half plane since (\begin{smallmatrix} 0 & 1 \\  1 & 0 \end{smallmatrix}) squares to identity. And around \pm \pi/2 something totally different (return to this later). The problem is now reduced to solving a non-homogeneous first order matrix equation of the form

\begin{aligned}\mathbf{z}' = \Omega \mathbf{z} + \mathbf{b}\end{aligned} \quad\quad\quad(136)

But we have the good fortune of being able to easily exponentiate and invert this matrix \Omega. The homogeneous problem

\begin{aligned}\mathbf{z}' = \Omega \mathbf{z}\end{aligned} \quad\quad\quad(137)

has the solution

\begin{aligned}\mathbf{z}_h(t) = e^{\Omega t} \mathbf{z}_{t=0}\end{aligned} \quad\quad\quad(138)

Assuming a specific solution z = e^{\Omega t}f(t) for the non-homogeneous equation, one finds z = \Omega^{-1} (e^{\Omega t} - I) \mathbf{b}. The complete solution with both the homogeneous and non-homogeneous parts is thus

\begin{aligned}\mathbf{z}(t) = e^{\Omega t} \mathbf{z}_0 + \Omega^{-1} (e^{\Omega t} - I) \mathbf{b}\end{aligned} \quad\quad\quad(139)

Going back to the pendulum problem, lets write

\begin{aligned}\omega = \sqrt{\frac{g}{l} {\left\lvert{\cos\theta_0}\right\rvert}}\end{aligned} \quad\quad\quad(140)

Below the horizontal we have

\begin{aligned}\Omega&= \omega\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix} \\ \Omega^{-1}&= -\frac{1}{{\omega}} \begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix} \\ e^{\Omega t}&=\cos(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sin(\omega t)\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\end{aligned} \quad\quad\quad(141)

Whereas above the horizontal we have

\begin{aligned}\Omega&= \omega\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix} \\ \Omega^{-1}&= \frac{1}{{\omega}} \begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix} \\ e^{\Omega t}&=\cosh(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sinh(\omega t)\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}\end{aligned} \quad\quad\quad(144)

In both cases we have

\begin{aligned}\begin{bmatrix}p_\theta \\ \alpha\end{bmatrix}&=\begin{bmatrix}m l^2 \omega & 0 \\ 0 & 1\end{bmatrix}\mathbf{z} \\ \mathbf{b} &= - \frac{g}{l}\frac{\sin\theta_0}{\omega}\begin{bmatrix}1 \\ 0\end{bmatrix}\end{aligned} \quad\quad\quad(147)

(where the real angle was \theta = \theta_0 + \alpha). Since in this case \Omega^{-1} and e^{\Omega t} commute, we have below the horizontal

\begin{aligned}\mathbf{z}(t)&= e^{\Omega t} (\mathbf{z}_0 - \Omega^{-1} \mathbf{b}) - \Omega^{-1} \mathbf{b} \\ &=\left(\cos(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sin(\omega t)\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\right)\left(\mathbf{z}_0 +\frac{1}{{\omega}} \begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\mathbf{b} \right)+\frac{1}{{\omega}} \begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\mathbf{b} \end{aligned}

Expanding out the \mathbf{b} terms and doing the same for above the horizontal we have respectively (below and above)

\begin{aligned}\mathbf{z}_\text{low}(t)&=\left(\cos(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sin(\omega t)\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\right)\left(\mathbf{z}_0 -\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\right)-\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix} \\ \mathbf{z}_\text{high}(t)&=\left(\cosh(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sinh(\omega t)\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}\right)\left(\mathbf{z}_0 +\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\right)+\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\end{aligned} \quad\quad\quad(149)

The only thing that is really left is re-insertion of the original momentum and position variables using the inverse relation

\begin{aligned}\mathbf{z} &=\begin{bmatrix}1/(m l^2 \omega) & 0 \\ 0 & 1\end{bmatrix}\begin{bmatrix}p_\theta \\ \theta - \theta_0\end{bmatrix}\end{aligned} \quad\quad\quad(151)

Will that final insertion do anything more than make things messier? Observe that the \mathbf{z}_0 only has a momentum component when expressed back in terms of the total angle \theta. Also recall that p_\theta = m l^2 \dot{\theta}, so we have

\begin{aligned}\mathbf{z} &=\begin{bmatrix}\dot{\theta}/\omega \\ \theta - \theta_0\end{bmatrix} \\ \mathbf{z}_0&=\begin{bmatrix}\dot{\theta}_{t=0}/\omega \\ 0\end{bmatrix} \\ \end{aligned} \quad\quad\quad(152)

If this is somehow mystically free of all math mistakes then we have the final solution

\begin{aligned}{\begin{bmatrix}\dot{\theta}(t)/\omega \\ \theta(t) - \theta_0\end{bmatrix}}_\text{low}&=\left(\cos(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sin(\omega t)\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix}\right)\left(\frac{\dot{\theta}(0)}{\omega}\begin{bmatrix}1 \\ 0\end{bmatrix}-\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\right)-\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix} \\ {\begin{bmatrix}\dot{\theta}(t)/\omega \\ \theta(t) - \theta_0\end{bmatrix}}_\text{high}&=\left(\cosh(\omega t)\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}+\sinh(\omega t)\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}\right)\left(\frac{\dot{\theta}(0)}{\omega}\begin{bmatrix}1 \\ 0\end{bmatrix}+\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\right)+\frac{g}{l}\frac{\sin\theta_0}{\omega^2} \begin{bmatrix}0 \\ 1 \end{bmatrix}\end{aligned} \quad\quad\quad(155)

A qualification is required to call this a solution since it is only a solution is the restricted range where \theta is close enough to \theta_0 (in some imprecisely specified sense). One could conceivably apply this in a recursive fashion however, calculating the result for a small incremental change, yielding the new phase space point, and repeating at the new angle.

The question of what the form of the solution in the neighborhood of \pm \pi/2 has also been ignored. That’s probably also worth considering but I don’t feel like trying now.

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: