Two Body Problem¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
Esercizio del libro¶
Introduzione¶
Consideriamo il problema a due corpi (ad esempio Terra–Luna) soggetti alla gravità Newtoniana. Indichiamo le masse come $m_1$ e $m_2$, le posizioni $\mathbf{r}_1, \mathbf{r}_2$ e i momenti $\mathbf{p}_1, \mathbf{p}_2$. La costante gravitazionale è $G$.
Hamiltoniana a due corpi¶
La Hamiltoniana totale è $$\large H(\mathbf{r}_1,\mathbf{r}_2,\mathbf{p}_1,\mathbf{p}_2) = \frac{\mathbf{p}_1^2}{2 m_1} + \frac{\mathbf{p}_2^2}{2 m_2} - \frac{G m_1 m_2}{|\mathbf{r}_1 - \mathbf{r}_2|}. $$
Trasformazione in coordinate centro di massa e relative¶
Definiamo il centro di massa e le coordinate relative: \begin{align} \mathbf{R} &= \frac{m_1 \mathbf{r}_1 + m_2 \mathbf{r}_2}{m_1 + m_2}, & \mathbf{r} &= \mathbf{r}_1 - \mathbf{r}_2,\\ \mathbf{P} &= \mathbf{p}_1 + \mathbf{p}_2, & \mathbf{p} &= \frac{m_2 \mathbf{p}_1 - m_1 \mathbf{p}_2}{m_1 + m_2}. \end{align}
Nel sistema centro di massa ($\mathbf{R}=0, \mathbf{P}=0$), la Hamiltoniana si riduce a: $$\large H(\mathbf{r},\mathbf{p}) = \frac{\mathbf{p}^2}{2\mu} - \frac{G m_1 m_2}{|\mathbf{r}|}, $$ dove $\mu = \frac{m_1 m_2}{m_1+m_2}$ è la massa ridotta.
Leggi di conservazione¶
L’energia e il momento angolare relativi sono: \begin{align} E &= \frac{\mathbf{p}^2}{2\mu} - \frac{G m_1 m_2}{r}, \\ \mathbf{L} &= \mathbf{r} \times \mathbf{p}. \end{align}
Per la Terra–Luna, valori tipici: $$\large E \sim -4\times 10^{28}\ \mathrm{J}, \quad |\mathbf{L}| \sim 3\times 10^{34}\ \mathrm{kg\,m^2/s}. $$
Hamilton–Jacobi separabile (problema di Keplero)¶
La funzione principale di Hamilton–Jacobi $S(\mathbf{r},t)$ soddisfa: $$\large \frac{\partial S}{\partial t} + H\!\left(\mathbf{r}, \nabla S \right) = 0. $$
Per il problema stazionario separabile: $$\large S(\mathbf{r},t) = W(\mathbf{r}) - E t, \quad \frac{1}{2\mu} (\nabla W)^2 - \frac{G m_1 m_2}{r} = E. $$
In coordinate polari $(r,\phi)$, con $L$ costante del moto: $$\large W(r,\phi) = W_r(r) + L \phi, \quad \left( \frac{d W_r}{d r} \right)^2 = 2\mu \left(E + \frac{G m_1 m_2}{r} \right) - \frac{L^2}{r^2}. $$
La soluzione formale radiale è: $$\large W_r(r) = \int_{r_0}^{r} \sqrt{ 2 \mu \left(E + \frac{G m_1 m_2}{r'} \right) - \frac{L^2}{r'^2} } \, dr'. $$
Quindi la funzione principale Hamilton–Jacobi completa è: $$\large \boxed{ S(r,\phi,t) = W_r(r) + L \phi - E t. } $$
Conclusioni¶
- Il problema a due corpi può essere ridotto a una massa ridotta $\mu$ in un potenziale centrale.
- Energia e momento angolare sono conservati.
- La separazione di Hamilton–Jacobi in coordinate polari porta direttamente alla soluzione Keplero ellittica.
- Questa formulazione permette di calcolare la funzione principale $S(r,\phi,t)$, base per le variabili di azione–angolo.
import numpy as np
from scipy.integrate import solve_ivp
def two_body_hamiltonian(p1, p2, r1, r2, m1, m2, G):
'''
Calculate Hamiltonian for the tujo-body problem.
:param p1: Momentum of first body.
:param p2: Momentum of second body.
:param r1: Position of first body.
:param r2: Position of second body.
:param m1: Mass of first body.
:param m2: Mass of second body.
:param G: Gravitational constant.
:return:
Hamiltonian value.
'''
T1 = np.dot(p1, p1) / (2 * m1)
T2 = np.dot(p2, p2) / (2 * m2)
V = -G * m1 * m2 / np.linalg.norm(r1 - r2)
return T1 + T2 + V
def transform_to_com_and_relative(r1, r2, p1, p2, m1, m2):
'''
Transform coordinates to center of mass and relative
coordinates.
:param r1: Position of first body.
:param r2: Position of second body.
:param p1: Momentum of first body.
:param p2: Momentum of second body.
:param m1: Mass of first body.
:param m2: Mass of second body.
:return:
Transformed positions and momenta.
'''
R = (m1 * r1 + m2 * r2) / (m1 + m2)
r = r1 - r2
P = p1 + p2
p = (m2 * p1 - m1 * p2) / (m1 + m2)
return R, r, P, p
def conservation_laws(r, p, ml, m2, G):
'''
Compute conservation laws: energy and angular momentum.
:param r: Relative position vector.
:param p: Relative momentum vector.
:param m1: Mass of first body.
:param m2: Mass of second body.
:param G: Gravitational constant.
:return:
Energy and angular momentum values.
'''
mu = m1 * m2 / (m1 + m2)
E = np.dot(p, p) / (2 * mu) - G * m1 * m2 / np.linalg.norm(r)
L = np.cross(r, p)
return E, L
def hamilton_jacobi_solution(r, p, m1, m2, G):
'''
Placeholder for Hamilton-Jacobi solution.
:param r: Relative position vector.
:param p: Relative momentum vector.
:param m1: Mass of first body.
:param m2: Mass of second body.
:return:
Principal function value (dummy implementation).
'''
# This is a simplified placeholder function
# Calculations for principal function go here
S = np.linalg.norm(p) * np.linalg.norm(r)
return S
def kepler_hj_radial_integrand(r, E, L, mu, k):
val = 2 * mu * (E + k / r) - L**2 / r**2
return np.sqrt(np.maximum(val, 0.0))
from scipy.integrate import quad
def kepler_Wr(r, E, L, mu, k, r0=1e7):
r = float(r) # forza scalare
Wr, _ = quad(kepler_hj_radial_integrand, r0, r, args=(E, L, mu, k))
return Wr
def hamilton_jacobi_kepler(r, phi, t, E, L, mu, k):
Wr = kepler_Wr(r, E, L, mu, k)
return Wr + L * phi - E * t
# Constants
G = 6.67430e-11 # Gravitational constant [m^3 kg^-1 s^-2]
m1 = 5.972e24 # Earth mass [kg]
m2 = 7.348e22 # Moon mass [kg]
mu = m1 * m2 / (m1 + m2)
k = G * m1 * m2
# Initial positions (Earth–Moon distance)
r1 = np.array([0.0, 0.0, 0.0])
r2 = np.array([384400000.0, 0.0, 0.0]) # [m]
# Orbital velocity of the Moon (approximately circular)
v2 = np.array([0.0, 1022.0, 0.0]) # [m/s]
# Barycentric velocities
v1 = - (m2 / m1) * v2
# Canonical momenta
p1 = m1 * v1 # [kg m/s]
p2 = m2 * v2 # [kg m/s]
# Hamiltonian and transformations
H = two_body_hamiltonian(p1, p2, r1, r2, m1, m2, G)
R, r, P, p = transform_to_com_and_relative(r1, r2, p1, p2, m1, m2)
# Conservation laws
E, L_vec = conservation_laws(r, p, m1, m2, G)
L = np.linalg.norm(L_vec)
# Coordinate polari dal vettore relativo
r_norm = np.linalg.norm(r)
phi = np.arctan2(r[1], r[0]) # piano orbitale xy
t = 0.0 # tempo iniziale
# HamiIton-Jacobi method (simple placeholder)
S = hamilton_jacobi_solution(r, p, m1, m2, G)
print(f"Hamiltonian: {H:.2e}")
print(
f"Energy < 0 Elleptic Orbit: {E:.2e}, "
f"Angular Momentum: ({L:.2e})"
)
print(f"Hamilton-Jacobi Principal Function (Placeholder): {S:.2e}")
#-----------------------------------------------------------------------
S_kepler = hamilton_jacobi_kepler(
r=r_norm,
phi=phi,
t=t,
E=E,
L=L,
mu=mu,
k=k
)
print("")
print(
f"Energy < 0 Elliptic Orbit: {E:.2e}, "
f"|L|: {L:.2e}"
)
print(f"Hamilton–Jacobi Principal Function (Kepler): {S_kepler:.2e} J·s")
Hamiltonian: -3.73e+28 Energy < 0 Elleptic Orbit: -3.73e+28, Angular Momentum: (2.89e+34) Hamilton-Jacobi Principal Function (Placeholder): 2.89e+34 Energy < 0 Elliptic Orbit: -3.73e+28, |L|: 2.89e+34 Hamilton–Jacobi Principal Function (Kepler): 9.07e+34 J·s