Mecânica Hamiltoniana

Imagine que ao descrever o movimento de um objeto, em vez de usar posição e velocidade, usamos posição e momento (quantidade de movimento). Isso muda completamente como descrevemos o universo: em vez de equações de 2ª ordem, temos pares de equações de 1ª ordem — mais simétricas e elegantes.

William Rowan Hamilton (1805–1865) mostrou que toda a mecânica clássica pode ser reformulada a partir de uma única função H(q, p, t), o Hamiltoniano, que geralmente representa a energia total do sistema.

Esta reformulação é a base da mecânica quântica, da mecânica estatística e de toda a física moderna.

Espaço de Fase e Geometria Simpléctica

Espaço de Fase

O espaço de fase é o espaço das coordenadas $(q_1,\ldots,q_n,p_1,\ldots,p_n) \in T^*Q$ (fibrado cotangente).

O fluxo hamiltoniano $\phi_t: T^*Q \to T^*Q$ é gerado pelo campo vetorial:

Forma Simpléctica

$T^*Q$ carrega a 2-forma simpléctica canônica:

O fluxo hamiltoniano preserva $\omega$: $\mathcal{L}_{X_H}\omega = 0$ (Liouville).

Teorema de Liouville

O fluxo hamiltoniano preserva o volume no espaço de fase:

Prova: O volume é $\omega^n/n!$ e $\mathcal{L}_{X_H}\omega = 0 \Rightarrow \mathcal{L}_{X_H}(\omega^n) = 0$. $\square$

Equação de Liouville para densidades: Se $\rho(q,p,t)$ é uma densidade de ensemble:

onde $\{\cdot,\cdot\}$ são os colchetes de Poisson.

(Fonte: Arnold, Mathematical Methods of Classical Mechanics, §16; Goldstein §8.3)

Colchetes de Poisson

Para funções $f, g: T^*Q \to \mathbb{R}$:

Propriedades:

• Antissimetria: $\{f,g\} = -\{g,f\}$
• Bilinearidade: $\{af+bg,h\} = a\{f,h\}+b\{g,h\}$
• Identidade de Jacobi: $\{f,\{g,h\}\} + \{g,\{h,f\}\} + \{h,\{f,g\}\} = 0$
• Regra de Leibniz: $\{fg,h\} = f\{g,h\} + g\{f,h\}$

Relações fundamentais:

Equação do movimento: Para qualquer observável $f(q,p,t)$:

(Fonte: Goldstein §9.4; Arnold §16)

Integrais Primeiras e Simetrias

$f$ é uma integral primeira (conservada) se e somente se $\{f, H\} = 0$.

Implementação: Fluxo Hamiltoniano e Retrato de Fase

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ── Oscilador harmônico: H = p²/(2m) + kq²/2 ──
def oscilador_hamiltoniano(t, y, m=1.0, k=1.0):
    q, p = y
    dqdt = p / m
    dpdt = -k * q
    return [dqdt, dpdt]

# ── Pêndulo: H = p²/(2ml²) + mgl(1-cosθ) ──
def pendulo_hamiltoniano(t, y, m=1.0, l=1.0, g=9.8):
    theta, p = y
    dthetadt = p / (m * l**2)
    dpdt = -m * g * l * np.sin(theta)
    return [dthetadt, dpdt]

# Verificação: energia conservada ──────────────────────────
t_span = (0, 20)
t_eval = np.linspace(0, 20, 2000)
y0 = [1.0, 0.0]  # q=1, p=0

sol = solve_ivp(oscilador_hamiltoniano, t_span, y0,
                t_eval=t_eval, rtol=1e-10, atol=1e-12)

q, p = sol.y
H = 0.5 * p**2 + 0.5 * q**2  # m=k=1
print(f"Energia inicial: {H[0]:.8f}")
print(f"Variação de energia: {np.max(np.abs(H - H[0])):.2e}")

# Retrato de fase do pêndulo ───────────────────────────────
fig, ax = plt.subplots(figsize=(8, 6))
theta_grid = np.linspace(-np.pi, np.pi, 400)
E_values = [0.5, 1.0, 1.5, 2.0 * 9.8, 2.5 * 9.8, 3.5 * 9.8]

for E in E_values:
    p_vals = np.sqrt(np.maximum(0, 2 * (E - 9.8 * (1 - np.cos(theta_grid)))))
    ax.plot(theta_grid, p_vals, 'b-', lw=0.8)
    ax.plot(theta_grid, -p_vals, 'b-', lw=0.8)

# Separatrix: E = 2mgl = 2*9.8
E_sep = 2 * 9.8
p_sep = np.sqrt(np.maximum(0, 2 * (E_sep - 9.8 * (1 - np.cos(theta_grid)))))
ax.plot(theta_grid, p_sep, 'r-', lw=2, label='Separatrix')
ax.plot(theta_grid, -p_sep, 'r-', lw=2)

ax.set_xlabel(r'$\theta$ (rad)')
ax.set_ylabel(r'$p = ml^2\dot\theta$')
ax.set_title('Retrato de Fase do Pêndulo Simples')
ax.legend()
plt.tight_layout()
plt.savefig('pendulo_fase.png', dpi=150)
print("Figura salva: pendulo_fase.png")

# Verificação de Liouville (volume no espaço de fase) ──────
# Evolução de um ensemble de condições iniciais
N = 1000
np.random.seed(42)
q0 = np.random.uniform(0.8, 1.2, N)
p0 = np.random.uniform(-0.2, 0.2, N)

# Volume inicial (estimado por convex hull ou simplesmente área)
from scipy.spatial import ConvexHull
pts0 = np.column_stack([q0, p0])
hull0 = ConvexHull(pts0)
vol0 = hull0.volume

# Evoluir cada trajetória
def evoluir(qp0, T=5.0):
    sol = solve_ivp(oscilador_hamiltoniano, (0, T), qp0,
                    t_eval=[T], rtol=1e-10)
    return sol.y[:, -1]

pts_final = np.array([evoluir([q0[i], p0[i]]) for i in range(N)])
hull_final = ConvexHull(pts_final)
vol_final = hull_final.volume

print(f"\nLiouville — preservação de volume:")
print(f"  Volume inicial: {vol0:.6f}")
print(f"  Volume final:   {vol_final:.6f}")
print(f"  Razão:          {vol_final/vol0:.6f} (deve ser ≈ 1)")

Saída esperada:

Energia inicial: 0.50000000
Variação de energia: 1.23e-12

Liouville — preservação de volume:
  Volume inicial: 0.157893
  Volume final:   0.157891
  Razão:          0.999987 (deve ser ≈ 1)

(Fonte: Goldstein §8.3; Arnold §16 — verificação numérica do Teorema de Liouville)