Part II: Hamiltonian MechanicsTopic 1 of 3

The Legendre Transform

The mathematical operation that takes us from \(L(q, \dot{q})\) to \(H(q, p)\) — replacing velocities with momenta as independent variables. The same transform connects all thermodynamic potentials.

Historical Context

Adrien-Marie Legendre introduced this transform in 1787 while studying the geometry of curves. The idea is deceptively simple: instead of describing a convex function by its values, describe it by its tangent lines. This "dual description" appeared independently in thermodynamics (Gibbs, 1870s) and in mechanics (Hamilton, 1834). The mathematical theory was placed on firm ground by Fenchel and Rockafellar in the 20th century as part of convex analysis.

In mechanics, the Legendre transform is the bridge from the Lagrangian to the Hamiltonian formulation. It replaces the velocity \(\dot{q}\) (a tangent vector) with the momentum \(p\) (a cotangent vector), revealing the natural symplectic structure of phase space.

1. The Legendre Transform: Definition

Given a function \(f(x)\) that is convex (\(f''(x) > 0\)), its Legendre transform is a new function \(g(p)\) defined by:

\[g(p) = \sup_x \left[px - f(x)\right] = p\,x(p) - f(x(p))\]

where \(x(p)\) is determined by the condition \(p = f'(x)\) (i.e., the supremum is attained where the derivative equals \(p\)). The key idea: we are switching from describing \(f\) in terms of its argument \(x\) to describing it in terms of its slope \(p = f'(x)\).

The Involution Property

The Legendre transform is its own inverse. If \(g\) is the Legendre transform of \(f\), then \(f\) is the Legendre transform of \(g\):

\[f(x) = \sup_p \left[xp - g(p)\right]\]

This is the involution (or duality) property. It means no information is lost: the Legendre transform is a bijection between convex functions and their duals.

2. Geometric Interpretation

The Legendre transform has a beautiful geometric meaning. Consider the graph of\(y = f(x)\). At each point \(x\), draw the tangent line with slope \(p = f'(x)\). This line has equation:

\[y = px - g(p)\]

where \(g(p)\) is the \(y\)-intercept of the tangent line (with a sign flip). So\(g(p)\) encodes the same curve by its family of tangent lines rather than its points. For a convex function, these two descriptions are equivalent — the curve is the envelope of its tangent lines.

Why Convexity Matters

If \(f\) is not convex, the map \(x \mapsto p = f'(x)\) is not one-to-one: different \(x\) values can give the same slope \(p\). The Legendre transform is then not invertible. Convexity ensures \(f'' > 0\), so \(f'\) is strictly monotonic, and the transform is a well-defined bijection.

3. From \(L(q, \dot{q})\) to \(H(q, p)\)

In mechanics, we apply the Legendre transform with respect to the velocity variables\(\dot{q}_i\), holding \(q_i\) and \(t\) fixed. The conjugate momentum is:

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

The Hamiltonian is the Legendre transform of \(L\) with respect to \(\dot{q}\):

The Hamiltonian

\[H(q, p, t) = \sum_i p_i \dot{q}_i - L(q, \dot{q}, t)\]

where \(\dot{q}_i\) on the right side must be expressed in terms of \(q_i\) and\(p_i\) by inverting the relation \(p_i = \partial L/\partial \dot{q}_i\).

Convexity Condition in Mechanics

The Legendre transform is well-defined when the "Hessian matrix" is non-singular:

\[\det\left(\frac{\partial^2 L}{\partial \dot{q}_i \partial \dot{q}_j}\right) \neq 0\]

For standard mechanical systems where \(T = \frac{1}{2}\sum_{ij} M_{ij}(q)\dot{q}_i\dot{q}_j\), this Hessian is just the mass matrix \(M_{ij}\), which is positive definite. So the Legendre transform always works for conventional mechanics. It fails for some singular Lagrangians that arise in gauge theories (e.g., electrodynamics), requiring the Dirac-Bergmann constraint algorithm.

4. Worked Examples

Example 1: Free Particle

\(L = \frac{1}{2}m\dot{x}^2\)

\(p = \frac{\partial L}{\partial \dot{x}} = m\dot{x} \Rightarrow \dot{x} = p/m\)

\(H = p\dot{x} - L = p \cdot \frac{p}{m} - \frac{1}{2}m\left(\frac{p}{m}\right)^2 = \frac{p^2}{2m}\)

As expected: the Hamiltonian is the kinetic energy expressed in terms of momentum.

Example 2: Harmonic Oscillator

\(L = \frac{1}{2}m\dot{x}^2 - \frac{1}{2}kx^2\)

\(p = m\dot{x} \Rightarrow \dot{x} = p/m\)

\(H = \frac{p^2}{m} - \frac{1}{2}\frac{p^2}{m} + \frac{1}{2}kx^2 = \frac{p^2}{2m} + \frac{1}{2}kx^2\)

\(H = T + V\): the total energy, which is conserved since \(H\) has no explicit time dependence.

Example 3: Particle in Electromagnetic Field

\(L = \frac{1}{2}m|\dot{\mathbf{r}}|^2 + \frac{e}{c}\dot{\mathbf{r}} \cdot \mathbf{A} - e\phi\)

\(\mathbf{p} = m\dot{\mathbf{r}} + \frac{e}{c}\mathbf{A}\) (canonical momentum includes the vector potential!)

\(\dot{\mathbf{r}} = \frac{1}{m}\left(\mathbf{p} - \frac{e}{c}\mathbf{A}\right)\)

\(H = \frac{1}{2m}\left|\mathbf{p} - \frac{e}{c}\mathbf{A}\right|^2 + e\phi\)

The canonical momentum \(\mathbf{p}\) is not the mechanical momentum \(m\dot{\mathbf{r}}\). This distinction is crucial in quantum mechanics.

Example 4: Relativistic Particle

\(L = -mc^2\sqrt{1 - |\dot{\mathbf{r}}|^2/c^2}\)

\(\mathbf{p} = \frac{m\dot{\mathbf{r}}}{\sqrt{1 - v^2/c^2}} = \gamma m \dot{\mathbf{r}}\)

\(H = \sqrt{|\mathbf{p}|^2 c^2 + m^2 c^4}\)

The famous relativistic energy-momentum relation \(E^2 = p^2c^2 + m^2c^4\) emerges directly from the Legendre transform.

5. Thermodynamic Analogy

The Legendre transform plays exactly the same role in thermodynamics, connecting the different thermodynamic potentials. The internal energy \(U(S, V, N)\) is the fundamental potential, and all others are obtained by Legendre transforms:

Transform

Result

\(S \to T = \partial U/\partial S\)

\(F = U - TS\) (Helmholtz free energy)

\(V \to P = -\partial U/\partial V\)

\(H_{\text{th}} = U + PV\) (Enthalpy)

\(S \to T\) and \(V \to P\)

\(G = U - TS + PV\) (Gibbs free energy)

The analogy is exact:

Mechanics

\(L(q, \dot{q}) \to H(q, p)\)

\(p = \partial L / \partial \dot{q}\)

\(H = p\dot{q} - L\)

Thermodynamics

\(U(S, V) \to F(T, V)\)

\(T = \partial U / \partial S\)

\(F = U - TS\)

6. The Routhian: A Hybrid Approach

The Routhian (introduced by Edward Routh, 1877) is a partial Legendre transform: we transform only the cyclic coordinates (those that don't appear in \(L\)), leaving the others in Lagrangian form. For cyclic coordinates \(q_c\) with conjugate momenta \(p_c\):

\[R(q_{\text{nc}}, \dot{q}_{\text{nc}}, p_c) = L - \sum_{\text{cyclic}} p_c \dot{q}_c\]

The Routhian acts as a Lagrangian for the non-cyclic coordinates and as a (negative) Hamiltonian for the cyclic ones. This is especially useful for systems with symmetry, like the symmetric top or the Kepler problem.

Example: Central Force with Routhian

For the Kepler problem in polar coordinates, \(L = \frac{1}{2}m(\dot{r}^2 + r^2\dot{\theta}^2) - V(r)\). Since \(\theta\) is cyclic, \(p_\theta = mr^2\dot{\theta} = l\) (conserved).

The Routhian is:

\[R = L - p_\theta \dot{\theta} = \frac{1}{2}m\dot{r}^2 - \frac{l^2}{2mr^2} - V(r)\]

This is effectively a 1D Lagrangian for the radial motion, with \(l^2/(2mr^2)\)acting as the centrifugal barrier in the effective potential.

Python Simulation: Legendre Transform Visualization

We visualize the Legendre transform geometrically, demonstrate the involution property, compute the Hamiltonian for several systems, and show the Routhian reduction of the Kepler problem.

Legendre Transform: Geometry, Involution, and Mechanics

Python

script.py250 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.optimize import brentq

# ============================================================
# Part 1: Geometric visualization of the Legendre Transform
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# f(x) = x^2/2 (simplest convex function)
x = np.linspace(-3, 3, 500)
f = x**2 / 2

ax1 = axes[0]
ax1.plot(x, f, 'b-', linewidth=2.5, label='f(x) = x^2/2')

# Draw tangent lines for several slopes p
slopes = [-2, -1, 0, 1, 2]
colors = ['#ef4444', '#f59e0b', '#22c55e', '#3b82f6', '#8b5cf6']
for p_val, col in zip(slopes, colors):
    # At x_0 where f'(x_0) = p, we have x_0 = p
    x0 = p_val
    g_val = p_val * x0 - x0**2 / 2   # = p^2/2
    # Tangent: y = p*(x - x0) + f(x0) = px - g(p)
    y_tang = p_val * x - g_val
    ax1.plot(x, y_tang, '--', color=col, linewidth=1, alpha=0.7,
             label='p={}, g={:.1f}'.format(p_val, g_val))
    ax1.plot(x0, x0**2/2, 'o', color=col, markersize=8)

ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.set_title('Legendre Transform: Tangent Lines')
ax1.set_ylim(-2, 5)
ax1.legend(fontsize=7, loc='upper left')
ax1.grid(True, alpha=0.3)

# The Legendre transform g(p) = p^2/2
p_arr = np.linspace(-3, 3, 500)
g_arr = p_arr**2 / 2

ax2 = axes[1]
ax2.plot(p_arr, g_arr, 'r-', linewidth=2.5, label='g(p) = p^2/2')
ax2.plot(x, f, 'b--', linewidth=1.5, alpha=0.5, label='f(x) = x^2/2')
ax2.set_xlabel('p (slope)')
ax2.set_ylabel('g(p) (y-intercept)')
ax2.set_title('Legendre Transform of f(x) = x^2/2')
ax2.legend()
ax2.grid(True, alpha=0.3)

print("--- Example 1: f(x) = x^2/2 ---")
print("Legendre transform: g(p) = p^2/2")
print("Involution check: LT(g) should give back f")

# Non-trivial example: f(x) = e^x
ax3 = axes[2]
x2 = np.linspace(-2, 2, 500)
f2 = np.exp(x2)
ax3.plot(x2, f2, 'b-', linewidth=2.5, label='f(x) = e^x')

# LT: p = e^x => x = ln(p), g(p) = p*ln(p) - p
p2 = np.linspace(0.05, 7, 500)
g2 = p2 * np.log(p2) - p2
ax3.plot(p2, g2, 'r-', linewidth=2.5, label='g(p) = p*ln(p) - p')
ax3.set_xlabel('x or p')
ax3.set_ylabel('Function value')
ax3.set_title('Legendre Transform of f(x) = e^x')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_ylim(-3, 8)

print("\n--- Example 2: f(x) = e^x ---")
print("Legendre transform: g(p) = p*ln(p) - p  (for p > 0)")

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

# ============================================================
# Part 2: Hamiltonian for different systems
# ============================================================
from scipy.integrate import solve_ivp

fig2, axes2 = plt.subplots(2, 2, figsize=(14, 10))

# Harmonic oscillator: phase space
k_spring, m_ho = 2.0, 1.0
omega = np.sqrt(k_spring / m_ho)

ax4 = axes2[0, 0]
for E in [0.5, 1.0, 2.0, 4.0]:
    # H = p^2/(2m) + kx^2/2 = E => ellipse
    theta_arr = np.linspace(0, 2*np.pi, 500)
    x_ell = np.sqrt(2*E/k_spring) * np.cos(theta_arr)
    p_ell = np.sqrt(2*m_ho*E) * np.sin(theta_arr)
    ax4.plot(x_ell, p_ell, linewidth=1.5, label='E={:.1f}'.format(E))

ax4.set_xlabel('x')
ax4.set_ylabel('p')
ax4.set_title('Harmonic Oscillator Phase Space (H = const)')
ax4.set_aspect('equal')
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3)

# Pendulum: phase space with Hamiltonian contours
ax5 = axes2[0, 1]
g_val = 9.81
l_val = 1.0
theta_grid = np.linspace(-np.pi, np.pi, 300)
p_grid = np.linspace(-15, 15, 300)
TH, PP = np.meshgrid(theta_grid, p_grid)
# H = p^2/(2ml^2) - mgl*cos(theta)
H_grid = PP**2 / (2 * m_ho * l_val**2) - m_ho * g_val * l_val * np.cos(TH)

levels = np.linspace(-10, 30, 25)
ax5.contour(TH, PP, H_grid, levels=levels, colors='royalblue', linewidths=0.8)
# Separatrix: H = mgl
sep_E = m_ho * g_val * l_val
p_sep = np.sqrt(2 * m_ho * l_val**2 * (sep_E + m_ho * g_val * l_val * np.cos(theta_grid)))
ax5.plot(theta_grid, p_sep, 'r-', linewidth=2, label='Separatrix')
ax5.plot(theta_grid, -p_sep, 'r-', linewidth=2)
ax5.set_xlabel('theta')
ax5.set_ylabel('p_theta')
ax5.set_title('Pendulum Phase Space (H contours)')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Routhian: Kepler effective potential
ax6 = axes2[1, 0]
r_arr = np.linspace(0.3, 5, 500)
k_grav = 1.0

for l_ang in [0.3, 0.5, 0.8, 1.0, 1.5]:
    V_eff = -k_grav / r_arr + l_ang**2 / (2 * r_arr**2)
    ax6.plot(r_arr, V_eff, linewidth=1.5, label='l={:.1f}'.format(l_ang))

ax6.axhline(y=0, color='gray', linestyle=':', alpha=0.5)
ax6.set_xlabel('r')
ax6.set_ylabel('V_eff(r)')
ax6.set_title('Routhian: Effective Potential for Kepler')
ax6.set_ylim(-3, 2)
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)

# Verify involution numerically for f(x) = x^4/4 + x^2/2
ax7 = axes2[1, 1]
x3 = np.linspace(-2, 2, 200)
f3 = x3**4 / 4 + x3**2 / 2
# p = f'(x) = x^3 + x
p_of_x = x3**3 + x3
# g(p) = p*x - f(x)
g_of_p = p_of_x * x3 - f3

ax7.plot(x3, f3, 'b-', linewidth=2.5, label='f(x) = x^4/4 + x^2/2')
ax7.plot(p_of_x, g_of_p, 'r-', linewidth=2.5, label='g(p) = LT[f](p)')

# Inverse transform: should recover f
# For each x, compute LT[g](x) = sup_p [xp - g(p)]
# We parametrically have g(p(x)), so LT[g](x) = x*p(x) - g(p(x))
f_recovered = x3 * p_of_x - g_of_p
ax7.plot(x3, f_recovered, 'g--', linewidth=1.5, label='LT[g](x) (should = f)')

max_err = np.max(np.abs(f_recovered - f3))
ax7.set_xlabel('x or p')
ax7.set_ylabel('Function value')
ax7.set_title('Involution: LT(LT(f)) = f  (err={:.1e})'.format(max_err))
ax7.legend(fontsize=8)
ax7.grid(True, alpha=0.3)

print("\n--- Involution Check ---")
print("f(x) = x^4/4 + x^2/2")
print("Max error |LT(LT(f)) - f| = {:.2e}".format(max_err))

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

# ============================================================
# Part 3: Kepler orbit via Routhian reduction
# ============================================================
print("\n--- Routhian Reduction of Kepler Problem ---")

def kepler_radial(t, y, l_val):
    r, pr = y
    dr = pr / m_ho
    dpr = l_val**2 / (m_ho * r**3) - k_grav / r**2
    return [dr, dpr]

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

# Radial phase portrait
ax8 = axes3[0]
l_fixed = 0.8
for E in [-0.4, -0.3, -0.2, -0.1, 0.0]:
    # Find turning points where pr = 0
    V_eff_func = lambda r: -k_grav / r + l_fixed**2 / (2 * m_ho * r**2)
    # Solve radial motion
    r_min_guess = l_fixed**2 / (2 * k_grav)
    try:
        r_min = brentq(lambda r: V_eff_func(r) - E, 0.1, 5.0)
        r_max = brentq(lambda r: V_eff_func(r) - E, r_min + 0.01, 20.0)
        r_plot = np.linspace(r_min, r_max, 300)
        pr_plot = np.sqrt(np.maximum(2 * m_ho * (E - V_eff_func(r_plot)), 0))
        ax8.plot(r_plot, pr_plot, 'b-', linewidth=1.2)
        ax8.plot(r_plot, -pr_plot, 'b-', linewidth=1.2)
    except:
        pass

ax8.set_xlabel('r')
ax8.set_ylabel('p_r')
ax8.set_title('Radial Phase Portrait (l={:.1f})'.format(l_fixed))
ax8.grid(True, alpha=0.3)

# Full orbit reconstructed from Routhian
ax9 = axes3[1]
r0_k, pr0_k = 1.5, 0
E_k = pr0_k**2/(2*m_ho) - k_grav/r0_k + l_fixed**2/(2*m_ho*r0_k**2)

sol_k = solve_ivp(kepler_radial, [0, 50], [r0_k, pr0_k],
                   args=(l_fixed,), max_step=0.01, rtol=1e-12, atol=1e-14)

r_k = sol_k.y[0]
# Reconstruct theta: d(theta)/dt = l/(m*r^2)
theta_k = np.cumsum(l_fixed / (m_ho * r_k**2) * np.gradient(sol_k.t))

x_k = r_k * np.cos(theta_k)
y_k = r_k * np.sin(theta_k)
ax9.plot(x_k, y_k, 'b-', linewidth=0.8)
ax9.plot(0, 0, 'yo', markersize=10, label='Center')
ax9.set_xlabel('x')
ax9.set_ylabel('y')
ax9.set_title('Orbit Reconstructed from Routhian')
ax9.set_aspect('equal')
ax9.legend()
ax9.grid(True, alpha=0.3)

print("Radial energy: E = {:.4f}".format(E_k))
print("Angular momentum: l = {:.4f}".format(l_fixed))

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

print("\nSimulation complete. Three plots generated.")
print("Plot 1: Geometric Legendre transform and tangent lines")
print("Plot 2: Phase spaces, effective potential, involution check")
print("Plot 3: Kepler problem via Routhian reduction")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Summary

The Legendre transform replaces a function's argument with its derivative as the independent variable.
Geometrically, it re-encodes a curve via its tangent lines instead of its points.
Convexity ensures the transform is well-defined and invertible (involution property).
In mechanics: \(H(q,p) = p\dot{q} - L(q,\dot{q})\) where \(p = \partial L / \partial \dot{q}\).
The thermodynamic potentials (Helmholtz, Gibbs, enthalpy) are all Legendre transforms of internal energy.
The Routhian is a partial Legendre transform, useful for systems with cyclic coordinates.

← Part II Overview Next: Hamilton's Equations →