Hamilton-Jacobi Equations¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
Calculate Hamilton-Jacobi Equation¶
Potenziale armonico¶
Consideriamo una particella di massa $m$ in una dimensione soggetta al potenziale armonico $$\large V(q)=\frac{1}{2}m\omega^2 q^2 . $$
L'Hamiltoniana è $$\large H(q,p)=\frac{p^2}{2m}+\frac{1}{2}m\omega^2 q^2 . $$
Equazione di Hamilton--Jacobi¶
L'equazione di Hamilton--Jacobi è $$\large H\!\left(q,\frac{\partial S}{\partial q},t\right) +\frac{\partial S}{\partial t}=0 , $$ ovvero $$\large \frac{1}{2m}\left(\frac{\partial S}{\partial q}\right)^2 +\frac{1}{2}m\omega^2 q^2 +\frac{\partial S}{\partial t}=0 . $$
Separazione delle variabili¶
Poiché l'Hamiltoniana non dipende esplicitamente dal tempo, si cerca una soluzione del tipo $$\large S(q,t)=W(q)-Et . $$
Sostituendo: $$\large \frac{1}{2m}\left(\frac{dW}{dq}\right)^2 +\frac{1}{2}m\omega^2 q^2 = E . $$
Questa è l'equazione di Hamilton--Jacobi stazionaria.
Soluzione per $W(q)$}¶
Si ottiene $$\large \frac{dW}{dq} =\sqrt{2m\left(E-\frac{1}{2}m\omega^2 q^2\right)} . $$
Integrando: $$\large W(q)=\int \sqrt{2mE-m^2\omega^2 q^2}\,dq . $$
L'integrale fornisce $$\large W(q)=\frac{m\omega}{2} \left[ q\sqrt{\frac{2E}{m\omega^2}-q^2} +\frac{2E}{m\omega^2} \arcsin\!\left( \frac{m\omega q}{\sqrt{2mE}} \right) \right]. $$
Ricostruzione della dinamica¶
Il momento canonico è dato da $$\large p=\frac{\partial S}{\partial q} =\frac{dW}{dq} =\sqrt{2mE-m^2\omega^2 q^2}. $$
Dalla relazione canonica $$\large \frac{\partial S}{\partial E}=\text{costante} $$ si ottiene $$\large \frac{\partial S}{\partial E} =\frac{\partial W}{\partial E}-t =\beta , $$ dove $\beta$ è una costante.
Calcoliamo $$\large \frac{\partial W}{\partial E} =\int \frac{m}{\sqrt{2mE-m^2\omega^2 q^2}}\,dq . $$
Segue: $$\large t+\beta =\int \frac{m\,dq}{\sqrt{2mE-m^2\omega^2 q^2}} . $$
Ponendo $$\large q = A\sin\theta, \qquad A=\sqrt{\frac{2E}{m\omega^2}}, $$ si ha $$\large dq = A\cos\theta\, d\theta $$ e $$\large \sqrt{2mE-m^2\omega^2 q^2} = m\omega A\cos\theta . $$
L'integrale diventa $$\large t+\beta =\frac{1}{\omega}\int d\theta =\frac{\theta}{\omega}. $$
Pertanto $$\large \theta=\omega t+\phi , $$ con $\phi=\omega\beta$.
Sostituendo: $$\large q(t)=A\sin(\omega t+\phi), \qquad A=\sqrt{\frac{2E}{m\omega^2}} . $$
Cenno su Hamilton--Jacobi e caos¶
Se l'equazione di Hamilton--Jacobi ammette una soluzione completa separabile, il sistema è integrabile e non può essere caotico.
Il caos emerge quando:
- l'equazione di Hamilton--Jacobi non è separabile,
- non esiste una funzione generatrice globale,
- le variabili azione-angolo cessano di esistere globalmente.
In questo senso il caos è legato al fallimento globale dell'equazione di Hamilton--Jacobi.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# -----------------------------
# Parametri fisici
# -----------------------------
m = 1.0 # massa
omega = 1.0 # frequenza
E = 1.0 # energia totale
# Ampiezza dalla HJ:
# A = sqrt(2E / (m omega^2))
A = np.sqrt(2 * E / (m * omega**2))
# Fase iniziale
phi = 0.0
# -----------------------------
# Hamiltoniana
# H(q,p) = p^2/(2m) + 1/2 m omega^2 q^2
# -----------------------------
def hamiltonian(q, p):
return p**2 / (2 * m) + 0.5 * m * omega**2 * q**2
# -----------------------------
# Equazioni di Hamilton
# dq/dt = ∂H/∂p
# dp/dt = -∂H/∂q
# -----------------------------
def hamilton_eq(t, y):
q, p = y
dqdt = p / m
dpdt = -m * omega**2 * q
return [dqdt, dpdt]
# -----------------------------
# Condizioni iniziali
# -----------------------------
q0 = A * np.sin(phi)
p0 = m * omega * A * np.cos(phi)
y0 = [q0, p0]
# -----------------------------
# Integrazione numerica
# -----------------------------
t_span = (0.0, 10.0)
t_eval = np.linspace(*t_span, 2000)
sol = solve_ivp(
hamilton_eq,
t_span,
y0,
t_eval=t_eval,
rtol=1e-10,
atol=1e-12
)
q_num = sol.y[0]
p_num = sol.y[1]
# -----------------------------
# Soluzione analitica (HJ)
# -----------------------------
q_analytic = A * np.sin(omega * t_eval + phi)
p_analytic = m * omega * A * np.cos(omega * t_eval + phi)
# -----------------------------
# Grafici
# -----------------------------
plt.figure(figsize=(10, 3))
plt.subplot(1, 2, 1)
plt.plot(t_eval, q_num, label="q(t) numerico")
plt.plot(t_eval, q_analytic, "--", label="q(t) analitico (HJ)")
plt.xlabel("t")
plt.ylabel("q")
plt.legend()
plt.title("Dinamica")
plt.subplot(1, 2, 2)
plt.plot(q_num, p_num)
plt.xlabel("q")
plt.ylabel("p")
plt.title("Spazio delle fasi")
plt.tight_layout()
plt.show()
# -----------------------------
# Controllo conservazione energia
# -----------------------------
H_vals = hamiltonian(q_num, p_num)
print("Energia media:", np.mean(H_vals))
print("Deviazione max:", np.max(np.abs(H_vals - E)))
Energia media: 0.9999999998819649 Deviazione max: 2.3501289803107284e-10
from sympy import *
from sympy import init_printing
init_printing() # abilita output LaTeX in Jupyter
import sympy as sp
def hamilton_jacobi_equation(H, S, q, t):
"""
Symbolic Hamilton-Jacobi equation.
Parameters:
H : sympy expression H(q, p)
S : sympy expression S(q, t)
q : sympy symbol (coordinate)
t : sympy symbol (time)
Returns:
sympy expression representing: H(q, ∂S/∂q) + ∂S/∂t
"""
# p = ∂S/∂q
p_sub = sp.diff(S, q)
# Sostituzione p -> ∂S/∂q dentro H
H_sub = H.subs(sp.symbols('p'), p_sub)
# Derivata temporale di S
dS_dt = sp.diff(S, t)
# Equazione di Hamilton-Jacobi
HJ = sp.simplify(H_sub + dS_dt)
return HJ
# ============== ESEMPIO: OSCILLATORE ARMONICO =====================
# Definizione simbolica di variabili
q, t, m, omega = sp.symbols('q t m omega', positive=True)
p = sp.symbols('p')
# Hamiltoniano dell'oscillatore armonico
H = p**2/(2*m) + m*omega**2 * q**2 / 2
# Ansatz separabile S(q,t) = W(q) - E t
E = sp.symbols('E')
W = sp.Function('W')(q)
S = W - E*t
# Equazione HJ
HJ_eq = hamilton_jacobi_equation(H, S, q, t)
print("Equazione di Hamilton-Jacobi:")
#sp.pprint(HJ_eq)
#print("\nForma algebrica:")
#print(HJ_eq)
display(HJ_eq)
Equazione di Hamilton-Jacobi:
Solve Hamilton-Jacobi Equations¶
from sympy import *
from sympy import init_printing
init_printing() # abilita output LaTeX in Jupyter
import sympy as sp
# Definizione simboli
q, t = sp.symbols('q t')
m = 1
omega = 1
E = sp.symbols('E') # energia simbolica
def solve_hj_eq_separation(q, t, E, m=1, omega=1):
W = sp.Function('W')(q)
# Derivata di W rispetto a q
p = sp.diff(W, q)
# Equazione di Hamilton-Jacobi (per oscillatore armonico)
HJ_eq = sp.Eq(p**2 / (2*m) + 1/2 * m * omega**2 * q**2, E)
# Risolvi per dW/dq
dWdq_sol = sp.solve(HJ_eq, sp.diff(W, q))
# Integra per ottenere W(q)
W_sol = [sp.integrate(sol, q) for sol in dWdq_sol]
return W_sol
# Esempio d'uso
solutions = solve_hj_eq_separation(q, t, E)
display(solutions[0])
display(solutions[1])