Part II: Hamiltonian MechanicsTopic 2 of 3

Hamilton's Equations

The symmetric first-order equations that govern all of classical mechanics: \(\dot{q} = \partial H/\partial p\),\(\dot{p} = -\partial H/\partial q\). Together with Poisson brackets and canonical transformations, they reveal the algebraic heart of mechanics.

Historical Context

William Rowan Hamilton published his reformulation of mechanics in two landmark papers (1834-1835), "On a General Method in Dynamics." His key insight was to treat position and momentum as independent variables on equal footing, reducing the second-order Euler-Lagrange equations to a system of first-order equations. Carl Gustav Jacobi then developed the Hamilton-Jacobi equation (1837), showing how to solve Hamilton's equations by finding a single generating function.

Poisson had introduced his bracket in 1809, but its deep significance was recognized only after Hamilton's work. In the 20th century, Dirac showed that the Poisson bracket is the classical limit of the quantum commutator: \(\{f, g\} \to \frac{1}{i\hbar}[\hat{f}, \hat{g}]\). This correspondence is the foundation of canonical quantization.

1. Derivation of Hamilton's Equations

Starting from \(H(q, p, t) = \sum_i p_i \dot{q}_i - L(q, \dot{q}, t)\), we compute the total differential:

\[dH = \sum_i \dot{q}_i\,dp_i + \sum_i p_i\,d\dot{q}_i - \sum_i \frac{\partial L}{\partial q_i}dq_i - \sum_i \frac{\partial L}{\partial \dot{q}_i}d\dot{q}_i - \frac{\partial L}{\partial t}dt\]

Using \(p_i = \partial L/\partial \dot{q}_i\), the \(d\dot{q}_i\) terms cancel:

\[dH = \sum_i \dot{q}_i\,dp_i - \sum_i \frac{\partial L}{\partial q_i}dq_i - \frac{\partial L}{\partial t}dt\]

But since \(H = H(q, p, t)\), we also have:

\[dH = \sum_i \frac{\partial H}{\partial q_i}dq_i + \sum_i \frac{\partial H}{\partial p_i}dp_i + \frac{\partial H}{\partial t}dt\]

Comparing coefficients of \(dq_i\), \(dp_i\), and \(dt\):

Hamilton's Equations of Motion

\[\dot{q}_i = \frac{\partial H}{\partial p_i}\]

\[\dot{p}_i = -\frac{\partial H}{\partial q_i}\]

\[\frac{\partial H}{\partial t} = -\frac{\partial L}{\partial t}\]

These are \(2n\) first-order ODEs, equivalent to the \(n\) second-order Euler-Lagrange equations. The beautiful antisymmetry — a plus sign in the first equation, a minus sign in the second — is the hallmark of Hamiltonian mechanics.

2. Poisson Brackets

The Poisson bracket of two phase-space functions \(f(q, p, t)\) and\(g(q, p, t)\) is defined as:

\[\{f, g\} = \sum_{i=1}^n \left(\frac{\partial f}{\partial q_i}\frac{\partial g}{\partial p_i} - \frac{\partial f}{\partial p_i}\frac{\partial g}{\partial q_i}\right)\]

Fundamental Properties

Antisymmetry: \(\{f, g\} = -\{g, f\}\)
Linearity: \(\{af + bg, h\} = a\{f, h\} + b\{g, h\}\)
Leibniz rule: \(\{fg, h\} = f\{g, h\} + g\{f, h\}\)
Jacobi identity: \(\{f, \{g, h\}\} + \{g, \{h, f\}\} + \{h, \{f, g\}\} = 0\)

Fundamental Brackets

The canonical coordinates themselves satisfy:

\[\{q_i, q_j\} = 0, \quad \{p_i, p_j\} = 0, \quad \{q_i, p_j\} = \delta_{ij}\]

These are the classical analogs of the canonical commutation relations in quantum mechanics: \([\hat{q}_i, \hat{p}_j] = i\hbar\,\delta_{ij}\).

Equations of Motion in Bracket Form

Hamilton's equations can be written compactly as:

\[\dot{f} = \{f, H\} + \frac{\partial f}{\partial t}\]

for any phase-space function \(f\). In particular, \(\dot{q}_i = \{q_i, H\}\) and\(\dot{p}_i = \{p_i, H\}\). A quantity \(f\) is conserved if and only if\(\{f, H\} = 0\) and \(\partial f / \partial t = 0\).

3. Canonical Transformations

A canonical transformation is a change of phase-space coordinates \((q, p) \to (Q, P)\) that preserves the form of Hamilton's equations:

\[\dot{Q}_i = \frac{\partial K}{\partial P_i}, \quad \dot{P}_i = -\frac{\partial K}{\partial Q_i}\]

for some new Hamiltonian \(K(Q, P, t)\). Equivalently, a transformation is canonical if and only if it preserves all Poisson brackets:

\[\{Q_i, Q_j\} = 0, \quad \{P_i, P_j\} = 0, \quad \{Q_i, P_j\} = \delta_{ij}\]

Generating Functions

Canonical transformations can be generated by four types of generating functions:

Type 1: \(F_1(q, Q, t)\) with \(p = \partial F_1/\partial q\), \(P = -\partial F_1/\partial Q\)

Type 2: \(F_2(q, P, t)\) with \(p = \partial F_2/\partial q\), \(Q = \partial F_2/\partial P\)

Type 3: \(F_3(p, Q, t)\) with \(q = -\partial F_3/\partial p\), \(P = -\partial F_3/\partial Q\)

Type 4: \(F_4(p, P, t)\) with \(q = -\partial F_4/\partial p\), \(Q = \partial F_4/\partial P\)

The new Hamiltonian is \(K = H + \partial F/\partial t\). If \(F\) has no explicit time dependence, then \(K = H\) — the Hamiltonian is unchanged.

Example: Point Transformation

A simple coordinate change \(Q = Q(q, t)\) is canonical. The Type-2 generating function is \(F_2 = \sum_i P_i Q_i(q, t)\), which gives:

\(p_j = \sum_i P_i \frac{\partial Q_i}{\partial q_j}\) — the usual transformation law for conjugate momenta.

4. The Symplectic Structure

Hamilton's equations can be written in matrix form using the phase-space vector \(\mathbf{z} = (q_1, \ldots, q_n, p_1, \ldots, p_n)^T\):

\[\dot{\mathbf{z}} = J \cdot \nabla_z H, \quad J = \begin{pmatrix} 0 & I_n \\ -I_n & 0 \end{pmatrix}\]

The matrix \(J\) is the symplectic matrix. It satisfies\(J^2 = -I_{2n}\) and \(J^T = -J\). A transformation is canonical if and only if its Jacobian matrix \(M\) is symplectic: \(M^T J M = J\).

The symplectic structure gives phase space a natural geometry — symplectic geometry — that is fundamentally different from Riemannian geometry. The key object is the symplectic 2-form \(\omega = \sum_i dp_i \wedge dq_i\), which is closed (\(d\omega = 0\)) and non-degenerate. This structure is preserved by Hamiltonian evolution.

5. The Hamilton-Jacobi Equation

The Hamilton-Jacobi equation is the ultimate goal of Hamiltonian mechanics: find a canonical transformation that makes the new Hamiltonian \(K = 0\), so that all new coordinates and momenta are constants of motion.

We seek a Type-2 generating function \(S(q, P, t)\) (called Hamilton's principal function) such that \(K = H + \partial S/\partial t = 0\). Since \(p_i = \partial S/\partial q_i\), we get:

The Hamilton-Jacobi Equation

\[H\left(q_1, \ldots, q_n, \frac{\partial S}{\partial q_1}, \ldots, \frac{\partial S}{\partial q_n}, t\right) + \frac{\partial S}{\partial t} = 0\]

This is a first-order partial differential equation for \(S(q, t)\). A complete solution (containing \(n\) free constants \(\alpha_1, \ldots, \alpha_n\) besides the trivial additive constant) gives the complete solution to the mechanical problem.

Example: Free Particle

\(H = p^2/(2m)\). The HJ equation is \(\frac{1}{2m}\left(\frac{\partial S}{\partial x}\right)^2 + \frac{\partial S}{\partial t} = 0\).

Try \(S = \alpha x - \frac{\alpha^2}{2m}t\). Then \(p = \partial S/\partial x = \alpha\) and\(\frac{\partial S}{\partial \alpha} = x - \frac{\alpha}{m}t = \beta\) (constant).

This gives \(x = \beta + \frac{\alpha}{m}t\) — uniform motion with momentum \(\alpha\) and initial position \(\beta\).

Connection to Quantum Mechanics

The Hamilton-Jacobi equation is the classical limit of the Schrodinger equation. Writing \(\psi = A\,e^{iS/\hbar}\) and taking \(\hbar \to 0\), the Schrodinger equation reduces to the Hamilton-Jacobi equation. The wavefronts of \(\psi\)become the surfaces of constant \(S\), and the rays (particle trajectories) are perpendicular to these surfaces.

Python Simulation: Hamiltonian Dynamics and Phase Portraits

We solve Hamilton's equations for several systems, visualize phase portraits, verify Poisson bracket relations numerically, and demonstrate canonical transformations.

Hamilton Equations: Phase Portraits and Poisson Brackets

Python

script.py295 lines

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Part 1: Phase portraits for different Hamiltonians
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# (a) Harmonic oscillator: H = p^2/(2m) + kx^2/2
ax1 = axes[0, 0]
m_val, k_val = 1.0, 4.0

def ho_hamilton(t, z):
    q, p = z
    dq = p / m_val            # dH/dp
    dp = -k_val * q           # -dH/dq
    return [dq, dp]

for E in [0.5, 1.0, 2.0, 4.0, 8.0]:
    q0 = np.sqrt(2 * E / k_val)
    sol = solve_ivp(ho_hamilton, [0, 10], [q0, 0], max_step=0.01)
    ax1.plot(sol.y[0], sol.y[1], 'b-', linewidth=1)

ax1.set_xlabel('q')
ax1.set_ylabel('p')
ax1.set_title('Harmonic Oscillator Phase Portrait')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# (b) Double well: H = p^2/2 - q^2/2 + q^4/4
ax2 = axes[0, 1]

def double_well(t, z):
    q, p = z
    dq = p
    dp = q - q**3
    return [dq, dp]

q_grid = np.linspace(-2, 2, 300)
p_grid = np.linspace(-2, 2, 300)
Q, P = np.meshgrid(q_grid, p_grid)
H_dw = P**2/2 - Q**2/2 + Q**4/4

levels = np.linspace(-0.3, 2, 30)
ax2.contour(Q, P, H_dw, levels=levels, colors='royalblue', linewidths=0.7)
# Separatrix at H = 0
ax2.contour(Q, P, H_dw, levels=[0], colors='red', linewidths=2)
ax2.set_xlabel('q')
ax2.set_ylabel('p')
ax2.set_title('Double Well: H = p^2/2 - q^2/2 + q^4/4')
ax2.grid(True, alpha=0.3)

# (c) Pendulum: H = p^2/(2ml^2) - mgl*cos(q)
ax3 = axes[1, 0]
m_p, l_p, g_p = 1.0, 1.0, 9.81

def pendulum_ham(t, z):
    q, p = z
    dq = p / (m_p * l_p**2)
    dp = -m_p * g_p * l_p * np.sin(q)
    return [dq, dp]

theta_grid = np.linspace(-np.pi, np.pi, 300)
pth_grid = np.linspace(-20, 20, 300)
TH, PTH = np.meshgrid(theta_grid, pth_grid)
H_pend = PTH**2 / (2*m_p*l_p**2) - m_p*g_p*l_p*np.cos(TH)

levels_p = np.linspace(-10, 40, 35)
ax3.contour(TH, PTH, H_pend, levels=levels_p, colors='royalblue', linewidths=0.7)
# Separatrix
H_sep = m_p * g_p * l_p
ax3.contour(TH, PTH, H_pend, levels=[H_sep], colors='red', linewidths=2)
ax3.set_xlabel('theta')
ax3.set_ylabel('p_theta')
ax3.set_title('Pendulum Phase Portrait (Hamiltonian)')
ax3.grid(True, alpha=0.3)

# (d) Henon-Heiles: H = (px^2 + py^2)/2 + (x^2 + y^2)/2 + x^2*y - y^3/3
ax4 = axes[1, 1]

def henon_heiles(t, z):
    x, y, px, py = z
    dxdt = px
    dydt = py
    dpx = -x - 2*x*y
    dpy = -y - x**2 + y**2
    return [dxdt, dydt, dpx, dpy]

# Plot Poincare section: y=0, py>0
x_sec, px_sec = [], []
for E_hh in [1.0/12]:
    for x0 in np.linspace(-0.4, 0.4, 30):
        for px0 in np.linspace(-0.4, 0.4, 30):
            # y=0, compute py from energy
            py2 = 2*E_hh - px0**2 - x0**2
            if py2 < 0:
                continue
            py0 = np.sqrt(py2)
            sol = solve_ivp(henon_heiles, [0, 200], [x0, 0, px0, py0],
                           max_step=0.05, rtol=1e-10, atol=1e-12)
            # Find y=0 crossings with py>0
            y_arr = sol.y[1]
            for j in range(1, len(y_arr)):
                if y_arr[j-1] < 0 and y_arr[j] >= 0:
                    # Linear interpolation
                    frac = -y_arr[j-1] / (y_arr[j] - y_arr[j-1])
                    x_cross = sol.y[0][j-1] + frac * (sol.y[0][j] - sol.y[0][j-1])
                    px_cross = sol.y[2][j-1] + frac * (sol.y[2][j] - sol.y[2][j-1])
                    x_sec.append(x_cross)
                    px_sec.append(px_cross)

ax4.plot(x_sec, px_sec, 'b.', markersize=0.3, alpha=0.5)
ax4.set_xlabel('x')
ax4.set_ylabel('px')
ax4.set_title('Henon-Heiles Poincare Section (E=1/12)')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()

print("--- Phase Portraits ---")
print("Generated 4 phase portraits:")
print("  (a) Harmonic oscillator (closed ellipses)")
print("  (b) Double well (separatrix + librations)")
print("  (c) Pendulum (librations + rotations)")
print("  (d) Henon-Heiles Poincare section (regular islands + chaos)")

# ============================================================
# Part 2: Poisson bracket verification
# ============================================================
print("\n--- Poisson Bracket Verification ---")

# Numerical Poisson bracket via finite differences
def poisson_bracket_num(f, g, q, p, eps=1e-7):
    """Compute {f, g} numerically at point (q, p)."""
    df_dq = (f(q + eps, p) - f(q - eps, p)) / (2 * eps)
    df_dp = (f(q, p + eps) - f(q, p - eps)) / (2 * eps)
    dg_dq = (g(q + eps, p) - g(q - eps, p)) / (2 * eps)
    dg_dp = (g(q, p + eps) - g(q, p - eps)) / (2 * eps)
    return df_dq * dg_dp - df_dp * dg_dq

# For harmonic oscillator H = p^2/(2m) + kq^2/2
H_ho = lambda q, p: p**2 / (2*m_val) + k_val * q**2 / 2
q_func = lambda q, p: q
p_func = lambda q, p: p
L_z = lambda q, p: q * p   # "angular momentum" analog

print("At (q, p) = (1.0, 2.0):")
print("  {{q, p}} = {:.10f}  (should be 1)".format(
    poisson_bracket_num(q_func, p_func, 1.0, 2.0)))
print("  {{q, H}} = {:.10f}  (should be p/m = {:.1f})".format(
    poisson_bracket_num(q_func, H_ho, 1.0, 2.0), 2.0/m_val))
print("  {{p, H}} = {:.10f}  (should be -kq = {:.1f})".format(
    poisson_bracket_num(p_func, H_ho, 1.0, 2.0), -k_val*1.0))

# ============================================================
# Part 3: Hamilton-Jacobi for harmonic oscillator
# ============================================================
fig2, axes2 = plt.subplots(1, 2, figsize=(14, 6))

# HJ solution: S(x, alpha, t) = alpha*x*cos(wt) - alpha^2*sin(2wt)/(4mw) ...
# For the HO, Hamilton's principal function separates as S = W(x) - Et
# W satisfies: (dW/dx)^2/(2m) + kx^2/2 = E
# => dW/dx = sqrt(2m(E - kx^2/2))
# => W = integral of that

omega_ho = np.sqrt(k_val / m_val)
ax5 = axes2[0]

for E_val in [1.0, 2.0, 4.0, 8.0]:
    x_max = np.sqrt(2 * E_val / k_val)
    x_range = np.linspace(-x_max * 0.999, x_max * 0.999, 500)
    dW = np.sqrt(np.maximum(2 * m_val * (E_val - k_val * x_range**2 / 2), 0))
    W = np.cumsum(dW) * (x_range[1] - x_range[0])
    ax5.plot(x_range, W, linewidth=1.5, label='E={:.0f}'.format(E_val))

ax5.set_xlabel('x')
ax5.set_ylabel("Hamilton's characteristic function W(x)")
ax5.set_title('Hamilton-Jacobi: W(x) for Harmonic Oscillator')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Canonical transformation: action-angle for HO
ax6 = axes2[1]
# In (q, p) space, HO trajectories are ellipses
# Canonical transform to (theta, J): q = sqrt(2J/(mw)) sin(theta), p = sqrt(2mwJ) cos(theta)
# Then H = omega * J (linear in J!)

t_aa = np.linspace(0, 10, 1000)
q0_aa, p0_aa = 2.0, 0.0
E_aa = k_val * q0_aa**2 / 2
J_aa = E_aa / omega_ho   # action variable

q_aa = q0_aa * np.cos(omega_ho * t_aa)
p_aa = -m_val * omega_ho * q0_aa * np.sin(omega_ho * t_aa)
theta_aa = np.mod(omega_ho * t_aa, 2 * np.pi)

ax6.plot(t_aa, theta_aa, 'b-', linewidth=1.5, label='theta(t) (angle variable)')
ax6.axhline(y=J_aa, color='r', linestyle='--', linewidth=2,
            label='J = E/omega = {:.2f} (action, const)'.format(J_aa))
ax6.set_xlabel('Time')
ax6.set_ylabel('Angle / Action')
ax6.set_title('Action-Angle Variables for HO')
ax6.legend()
ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output2.png', dpi=150, bbox_inches='tight')
plt.close()

# ============================================================
# Part 4: Symplectic integrator comparison
# ============================================================
print("\n--- Symplectic vs Non-Symplectic Integration ---")

fig3, axes3 = plt.subplots(1, 2, figsize=(14, 6))

# Kepler problem
def kepler_ham(t, z):
    x, y, px, py = z
    r = np.sqrt(x**2 + y**2)
    r3 = r**3
    return [px, py, -x/r3, -y/r3]

# Standard RK45 (non-symplectic)
T_final = 200
sol_rk = solve_ivp(kepler_ham, [0, T_final], [1, 0, 0, 0.8],
                    max_step=0.1, rtol=1e-8, atol=1e-10)

# Stormer-Verlet (symplectic)
dt_sv = 0.01
N_sv = int(T_final / dt_sv)
z_sv = np.zeros((4, N_sv + 1))
z_sv[:, 0] = [1, 0, 0, 0.8]
t_sv = np.linspace(0, T_final, N_sv + 1)

for i in range(N_sv):
    x, y, px, py = z_sv[:, i]
    # Half-step momentum
    r = np.sqrt(x**2 + y**2)
    r3 = r**3
    px_half = px - 0.5 * dt_sv * x / r3
    py_half = py - 0.5 * dt_sv * y / r3
    # Full-step position
    x_new = x + dt_sv * px_half
    y_new = y + dt_sv * py_half
    # Half-step momentum
    r_new = np.sqrt(x_new**2 + y_new**2)
    r3_new = r_new**3
    px_new = px_half - 0.5 * dt_sv * x_new / r3_new
    py_new = py_half - 0.5 * dt_sv * y_new / r3_new
    z_sv[:, i+1] = [x_new, y_new, px_new, py_new]

# Energy comparison
E_rk = 0.5*(sol_rk.y[2]**2 + sol_rk.y[3]**2) - 1.0/np.sqrt(sol_rk.y[0]**2 + sol_rk.y[1]**2)
E_sv = 0.5*(z_sv[2]**2 + z_sv[3]**2) - 1.0/np.sqrt(z_sv[0]**2 + z_sv[1]**2)

ax7 = axes3[0]
ax7.plot(sol_rk.t, E_rk - E_rk[0], 'r-', linewidth=0.8, label='RK45 (non-symplectic)', alpha=0.8)
ax7.plot(t_sv, E_sv - E_sv[0], 'b-', linewidth=0.8, label='Stormer-Verlet (symplectic)', alpha=0.8)
ax7.set_xlabel('Time')
ax7.set_ylabel('Energy drift (E - E0)')
ax7.set_title('Energy Conservation: Symplectic vs RK45')
ax7.legend()
ax7.grid(True, alpha=0.3)

# Orbit comparison
ax8 = axes3[1]
ax8.plot(z_sv[0, ::10], z_sv[1, ::10], 'b-', linewidth=0.5, alpha=0.7, label='Symplectic')
ax8.plot(sol_rk.y[0], sol_rk.y[1], 'r-', linewidth=0.5, alpha=0.7, label='RK45')
ax8.plot(0, 0, 'yo', markersize=8)
ax8.set_xlabel('x')
ax8.set_ylabel('y')
ax8.set_title('Kepler Orbit: Long-Time Behavior')
ax8.set_aspect('equal')
ax8.legend()
ax8.grid(True, alpha=0.3)

print("RK45 energy drift at t={}: {:.6e}".format(T_final, E_rk[-1] - E_rk[0]))
print("Verlet energy drift at t={}: {:.6e}".format(T_final, E_sv[-1] - E_sv[0]))

plt.tight_layout()
plt.savefig('output3.png', dpi=150, bbox_inches='tight')
plt.close()

print("\nSimulation complete. Three plots generated.")
print("Plot 1: Phase portraits (HO, double well, pendulum, Henon-Heiles)")
print("Plot 2: Hamilton-Jacobi function, action-angle variables")
print("Plot 3: Symplectic vs non-symplectic integration comparison")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Summary

Hamilton's equations are \(2n\) first-order ODEs treating \(q\) and \(p\) symmetrically.
Poisson brackets encode the algebraic structure of mechanics: \(\dot{f} = \{f, H\}\).
The fundamental brackets \(\{q_i, p_j\} = \delta_{ij}\) become commutators in quantum mechanics.
Canonical transformations preserve the Poisson bracket structure (symplecticity).
The Hamilton-Jacobi equation reduces mechanics to solving a single PDE for the generating function \(S\).
Symplectic integrators preserve phase-space structure and give superior long-time energy conservation.

← Legendre Transform Next: Phase Space & Liouville →