6.3. Liouville Theorem Visualization

HamiltonianPendulum class to generate basic pendulum plots, applied to problems from Taylor, but now solving using Hamilton’s equations.

• Last revised 16-Mar-2019 by Dick Furnstahl (furnstahl.1@osu.edu).

For a pendulum, the Lagrangian is

\(\begin{align} \mathcal{L} = \frac12 m L^2 \dot\phi^2 - mgL(1 - \cos\phi) \end{align}\)

which tells us that the generalized momentum is

\(\begin{align} \frac{\partial\mathcal{L}}{\partial \dot\phi} = m L^2 \dot\phi \equiv p_\phi \;. \end{align}\)

We can invert this equation to find \(\dot\phi = p_\phi / m L^2\). Constructing the Hamiltonian by Legendre transformation we find

\(\begin{align} \mathcal{H} &= \dot\phi p_\phi - \mathcal{L} \\ &= \frac{p_\phi^2}{m L^2} - \frac12 m L^2 \dot\phi^2 + mgL(1 - \cos\phi) \\ &= \frac{p_\phi^2}{2 m L^2} + mgL(1 - \cos\phi) \;. \end{align}\)

Thus \(\mathcal{H}\) is simply \(T + V\). Hamilton’s equations are

\(\begin{align} \dot\phi &= \frac{\partial\mathcal{H}}{\partial p_\phi} = \frac{p_\phi}{m L^2} \\ \dot p_\phi &= -\frac{\partial\mathcal{H}}{\partial \phi} = -mgL \sin\phi \;. \end{align}\)

A consequence of Liouville’s theorem is that if we consider a region in phase space as a set of initial conditions and then evolve each point in that space, the region can change shape but its volume will not change.

We will try this out for the simple pendulum by first picking an area (e.g., a rectangle) in \(p_\phi\) by \(\phi\) and then selecting random points inside of it. Each of the points is then propagated a time step by Hamilton’s equations and plotted again.

%matplotlib inline

import numpy as np
from numpy.random import uniform, rand

from scipy.integrate import solve_ivp
from scipy.spatial import ConvexHull

import matplotlib.pyplot as plt

# The dpi (dots-per-inch) setting will affect the resolution and how large
#  the plots appear on screen and printed.  So you may want/need to adjust 
#  the figsize when creating the figure.
plt.rcParams['figure.dpi'] = 100.    # this is the default for notebook

# Change the common font size (smaller when higher dpi)
font_size = 10
plt.rcParams.update({'font.size': font_size})

Pendulum class and utility functions

class HamiltonianPendulum():
    """
    Pendulum class implements the parameters and Hamilton's equations for 
     a simple pendulum (no driving or damping).
     
    Parameters
    ----------
    omega_0 : float
        natural frequency of the pendulum (\sqrt{g/l} where l is the 
        pendulum length) 
    mass : float
        mass of pendulum

    Methods
    -------
    dy_dt(t, y)
        Returns the right side of the differential equation in vector y, 
        given time t and the corresponding value of y.
    """
    def __init__(self, L=1., mass=1., g=1.
                ):
        self.L = L
        self.g = g
        self.omega_0 = np.sqrt(g/L)
        self.mass = mass
    
    def dy_dt(self, t, y):
        """
        This function returns the right-hand sides of Hamilton's eqs
        [dq/dt dp/dt] --> [dphi/dt dp_phi/dt]
        
        Parameters
        ----------
        t : float
            time 
        y : float
            A 2-component vector with y[0] = phi(t) and y[1] = p_phi(t)
            
        Returns
        -------
        
        """
        return [y[1] / (self.mass * self.L**2), 
                -self.mass * self.g * self.L * np.sin(y[0]) 
               ]
    
    
    def solve_ode(self, t_pts, q_0, p_0, 
                  abserr=1.0e-10, relerr=1.0e-10):
        """
        Solve the ODE given initial conditions.
        Specify smaller abserr and relerr to get more precision.
        """
        y = [q_0, p_0] 
        solution = solve_ivp(self.dy_dt, (t_pts[0], t_pts[-1]), 
                             y, t_eval=t_pts, 
                             atol=abserr, rtol=relerr)
        q, p = solution.y

        return q, p

def plot_y_vs_x(x, y, axis_labels=None, label=None, title=None, 
                color=None, linestyle=None, semilogy=False, loglog=False,
                ax=None):
    """
    Generic plotting function: return a figure axis with a plot of y vs. x,
    with line color and style, title, axis labels, and line label
    """
    if ax is None:        # if the axis object doesn't exist, make one
        ax = plt.gca()

    if (semilogy):
        line, = ax.semilogy(x, y, label=label, 
                            color=color, linestyle=linestyle)
    elif (loglog):
        line, = ax.loglog(x, y, label=label, 
                          color=color, linestyle=linestyle)
    else:
        line, = ax.plot(x, y, label=label, 
                    color=color, linestyle=linestyle)

    if label is not None:    # if a label if passed, show the legend
        ax.legend()
    if title is not None:    # set a title if one if passed
        ax.set_title(title)
    if axis_labels is not None:  # set x-axis and y-axis labels if passed  
        ax.set_xlabel(axis_labels[0])
        ax.set_ylabel(axis_labels[1])

    return ax, line

def start_stop_indices(t_pts, plot_start, plot_stop):
    start_index = (np.fabs(t_pts-plot_start)).argmin()  # index in t_pts array 
    stop_index = (np.fabs(t_pts-plot_stop)).argmin()  # index in t_pts array 
    return start_index, stop_index

Make Liouville theorem visualization

# Labels for individual plot axes
phi_vs_time_labels = (r'$t$', r'$\phi(t)$')
p_phi_vs_time_labels = (r'$t$', r'$p_\phi(t)$')
phase_space_labels = (r'$\phi$', r'$p_\phi$')

L = 1.
g = 1.
mass = 1.

# Instantiate a pendulum 
p1 = HamiltonianPendulum(L=L, g=g, mass=mass)

# Generate initial conditions
num_samples = 300
phi_and_p_low = [1.5, 0.0]
phi_and_p_high = [2.0, 0.5]
out = uniform(low=phi_and_p_low, high=phi_and_p_high, size=(num_samples,2))

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(out[:, 0], out[:, 1], '.')
ax.set_aspect(1)
ax.set_xlabel(r'$\phi$')
ax.set_ylabel(r'$p_\phi$')
ax.set_xlim(-3., 3.)
ax.set_ylim(-2., 2.)

fig.tight_layout()

fig = plt.figure()
ax = fig.add_subplot(1,1,1)

for phi, p in out:
    ax.plot(phi, p, '.')   

ax.set_xlabel(r'$\phi$')
ax.set_ylabel(r'$p_\phi$')
ax.set_xlim(-3., 3.)
ax.set_ylim(-2., 2.)
ax.set_aspect(1)

fig.tight_layout()

# Common plotting time (generate the full time then use slices)
delta_t = 0.5
t_start = 0.
num_t = 15
t_end = t_start + (num_t - 1) * delta_t  

t_pts = np.arange(t_start, t_end+delta_t, delta_t)  

# start the plot!
#fig = plt.figure(figsize=(15,5.*num_t/3.))
fig = plt.figure(figsize=(15,5.*num_t/3.))
overall_title = 'Liouville theorem visualization with pendulum'

plot_num = 0
all_orbits = []
for phi_0, p_phi_0 in out:
    phi, p_phi = p1.solve_ode(t_pts, phi_0, p_phi_0)
    all_orbits.append([phi, p_phi])
        
for i in np.arange(num_t):
    ax_c = fig.add_subplot(int(num_t/3), 3, i+1)
    ax_c.set_title(rf'$t = {t_pts[i]:.1f}$')
    for j in np.arange(num_samples):
        phi_now, p_phi_now = all_orbits[j]
        ax_c.plot(phi_now[i], p_phi_now[i], '.')
        ax_c.set_xlim(-3., 3.)
        ax_c.set_ylim(-2., 2.)    

fig.tight_layout()
#fig.suptitle(overall_title, va='top')

fig.savefig('Liouville_test.png', bbox_inches='tight')