Stability of Orbits¶

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

Esercizio del libro¶

In [29]:
import numpy as np

def hamiltonian(p_r, L, r, mu, V):
    '''
    Computes the Hamiltonian for a two-body central force problem.
        :param p_r: Radial momentum.
        :param L: Angular momentum.
        :param r: Radial distance.
        :param mu: Reduced mass.
        :param V: Potential energy function.
    :return: 
        Hamiltonian value.
    '''
    return (p_r**2 / (2 * mu)) + (L**2 / (2 * mu * r**2)) + V(r)

def perturb_orbit(r, theta, delta_r, delta_theta):
    '''
    Introduce perturbations to radial distance and angle.
        :param r: Radial distance.
        :param theta: Angular coordinate.
        :param delta_r: Perturbation in radial distance.
        :param delta_theta: Perturbation in angle.
    :return: 
        Perturbed values.
    '''
    return r + delta_r, theta + delta_theta

def runge_lenz_vector(p, L, mu, k, r):
    '''
    Calculate the Runge-Lenz vector A.
        :param p: Momentum vector.
        :param L: Angular momentum vector.
        :param mu: Reduced mass.
        :param k: Gravitational constant times mass product.
        :param r: Position vector.
    :return: 
        Runge-Lenz vector.
    '''
    return np.cross(p, L) - mu * k * (r / np.linalg.norm(r)) 

def stability_criteria(K):
    '''
    Assess stability based on the eigenvalues of the Hessian.
        :param K: Hessian matrix.
    :return: 
        Stability status (True if stable, False otherwise).
    '''
    eigenvalues = np.linalg.eigvals(K)
    return np.all(eigenvalues > 0)

def evaluate_stability_of_circular_orbit(r_0, p_0, L, mu, V, delta_r):
    '''
    Evaluates the stability of a circular orbit under perturbations.
        :param r_0: Initial radial distance.
        :param p_0: Initial momentum.
        :param L: Angular momentum.
        :param mu: Reduced mass.
        :param V: Potential energy function.
        :param delta_r: Small perturbation in radial distance.
    :return: 
        Stability status.
    '''
    H_0 = hamiltonian(p_0, L, r_0, mu, V)
    r_perturbed = r_0 + delta_r
    H_prime = hamiltonian(p_0, L, r_perturbed, mu, V)
    delta_H = H_prime - H_0
    print(delta_H)
    return delta_H > 0

# Example potential function for a Keplerian potential
def V_kepler(r):
    G = 6.674e-11 # gravitational constant
    M = 5.972e24 # mass of a celestial body
    return -G * M / r

# Define sample initial conditions
r_0 = 7000e3 # initial radial distance in meters
p_0 = 0 # initial momentum (assume circular orbit)
L = 1e10 # artificial angular momentum
mu = 1 # reduced mass for simplification
delta_r = 0.0001 # small perturbation 

# Calculate perturbations and stability 
r_prime, theta_prime = perturb_orbit(r_0, 0, delta_r, 0)

stable = evaluate_stability_of_circular_orbit(r_0, p_0, L, mu, V_kepler, delta_r)
print("Stability of perturbed orbit:", "Stable" if stable else "Unstable")

0.0007842555642127991
Stability of perturbed orbit: Stable

Stabilità di un'orbita circolare: chiarimenti concettuali¶

Consideriamo il moto di una particella di massa ridotta $m=1$ in un campo gravitazionale newtoniano a simmetria centrale.
L'Hamiltoniana ridotta è:

$$\large H = \frac{p_r^2}{2} + V_{\text{eff}}(r), \qquad V_{\text{eff}}(r) = -\frac{\mu}{r} + \frac{L^2}{2r^2}, $$

dove $L$ è il momento angolare, costante in assenza di coppie esterne.

Orbita circolare come punto fisso

Un'orbita circolare non è definita solo dal raggio, ma dalla terna:

$$\large (r_0,\; p_r=0,\; L=\sqrt{\mu r_0}), $$

che soddisfa:

$$\large \frac{dV_{\text{eff}}}{dr}(r_0)=0. $$

Essa rappresenta un $\textit{punto fisso}$ nello spazio delle fasi $(r,p_r)$.

Stabilità lineare (locale)

La stabilità dell'orbita circolare è determinata dal segno del secondo derivato del potenziale efficace:

$$\large V_{\text{eff}}''(r_0) = -\frac{2\mu}{r_0^3} + \frac{3L^2}{r_0^4} = \frac{\mu}{r_0^3} > 0. $$

Ne segue che il punto fisso è $\textit{ellittico}$ e quindi stabile nel senso di Lyapunov:
perturbazioni infinitesime producono oscillazioni radiali periodiche.

Significato fisico di una perturbazione $\delta r$

Imporre una perturbazione radiale

$$\large (r_0,\;0,\;L) \;\longrightarrow\; (r_0+\delta r,\;0,\;L) $$

non rappresenta una nuova orbita circolare, ma un $\textit{nuovo stato iniziale}$ ottenuto tramite un intervento esterno.

Fisicamente ciò equivale a:

  • applicare una forza radiale,
  • senza applicare coppie (quindi $L$ resta costante),
  • compiendo un lavoro esterno.

Lavoro esterno ed energia

Il lavoro necessario a imporre $\delta r$ è:

$$\large W(\delta r) = H(r_0+\delta r) - H(r_0) = V_{\text{eff}}(r_0+\delta r) - V_{\text{eff}}(r_0). $$

Il segno e l'entità di $W$ dipendono fortemente dalla direzione della perturbazione.

Asimmetria tra $\delta r>0$ e $\delta r<0$

Il potenziale efficace presenta una forte asimmetria:

  • per $\delta r>0$ il pozzo è largo e poco ripido;
  • per $\delta r<0$ la barriera centrifuga cresce rapidamente come $1/r^2$.

Ne consegue che:

  • grandi perturbazioni verso l'esterno restano legate;
  • perturbazioni relativamente piccole verso l'interno possono portare a caduta.

Questa non è un'instabilità lineare, ma una \emph{instabilità non lineare globale}.

Equazioni del moto e oscillazioni radiali

Le equazioni ridotte sono:

$$\large \dot r = p_r, \qquad \dot p_r = -\frac{dV_{\text{eff}}}{dr}. $$

Linearizzando attorno a $r_0$ si ottengono oscillazioni armoniche con frequenza:

$$\large \omega_r^2 = V_{\text{eff}}''(r_0), $$

che nel potenziale newtoniano coincide con la frequenza angolare: $\omega_r = \omega_\theta$.

Interpretazione fisica finale

  • L'orbita circolare è stabile solo $\textit{localmente}$.
  • Imporre $\delta r$ richiede sempre lavoro esterno.
  • La conservazione di $L$ è coerente con l'assenza di coppie.
  • L'instabilità verso l'interno riflette la struttura del potenziale efficace.

Questo stesso meccanismo spiega l'evoluzione secolare di satelliti come Phobos, la cui perdita di momento angolare porta alla scomparsa del punto fisso stesso.

Stabilità orbita circolare¶

In [39]:
import numpy as np
import matplotlib.pyplot as plt

# ------------------------
# Parametri fisici
# ------------------------
mu = 1.0      # GM
L  = 1.0      # momento angolare

# r orbita circolare
r0 = L**2 / mu

# ------------------------
# Potenziale efficace
# ------------------------
def V_eff(r):
    return -mu/r + L**2/(2*r**2)

def dV_dr(r):
    return mu/r**2 - L**2/r**3

def d2V_dr2(r):
    return -2*mu/r**3 + 3*L**2/r**4

# ------------------------
# Hessiano (1D)
# ------------------------
Hessian = d2V_dr2(r0)

print(f"r0 = {r0:.3f}")
print(f"Hessiano V_eff''(r0) = {Hessian:.3f}")

if Hessian > 0:
    print("Orbita circolare STABILE (lineare)")
else:
    print("Orbita circolare INSTABILE")

# ------------------------
# Test energetico con perturbazioni finite
# ------------------------
delta = np.linspace(-0.6*r0, 1.5*r0, 400)
r = r0 + delta

H0 = V_eff(r0)
H1 = V_eff(r)

# ------------------------
# Plot
# ------------------------
plt.figure()
plt.plot(r, H1, label="V_eff(r)")
plt.axhline(H0, linestyle="--", label="Energia orbita circolare")
plt.axvline(r0, linestyle=":", label="r0")
plt.xlabel("r")
plt.ylabel("Energia")
plt.legend()
plt.title("Stabilità orbitale: locale vs globale")
plt.show()

r0 = 1.000
Hessiano V_eff''(r0) = 1.000
Orbita circolare STABILE (lineare)
No description has been provided for this image

Oscillazioni Orbitali - Spazio delle fasi¶

In [61]:
import numpy as np
import matplotlib.pyplot as plt

# ------------------------
# Parametri fisici
# ------------------------
mu = 1.0
L  = 1.0

r0 = L**2 / mu          # r orbita circolare

# ------------------------
# Equazioni del moto
# ------------------------
def deriv(y):
    r, pr = y
    drdt  = pr
    dprdt = -(mu/r**2 - L**2/r**3)
    return np.array([drdt, dprdt])

# ------------------------
# Integrazione (RK4)
# ------------------------
def integrate(y0, dt, N):
    y = np.zeros((N, 2))
    y[0] = y0

    for i in range(N-1):
        k1 = deriv(y[i])
        k2 = deriv(y[i] + 0.5*dt*k1)
        k3 = deriv(y[i] + 0.5*dt*k2)
        k4 = deriv(y[i] + dt*k3)
        y[i+1] = y[i] + dt*(k1 + 2*k2 + 2*k3 + k4)/6

    return y

# ------------------------
# Condizioni iniziali
# ------------------------
delta_r = -0.45 * r0     # perturbazione piccola
y0 = np.array([r0 + delta_r, 0.0])

# ------------------------
# Time grid
# ------------------------
dt = 0.01
T  = 50
N  = int(T/dt)

sol = integrate(y0, dt, N)

t  = np.linspace(0, T, N)
r  = sol[:,0]
pr = sol[:,1]

# ------------------------
# Plot oscillazioni radiali
# ------------------------
plt.figure()
plt.plot(t, r, label="r(t)")
plt.axhline(r0, linestyle="--", label="r0")
plt.xlabel("tempo")
plt.ylabel("r")
plt.legend()
plt.title("Oscillazioni radiali attorno all'orbita circolare")
plt.show()

# ------------------------
# Ritratto di fase
# ------------------------
plt.figure()
plt.plot(r, pr)
plt.xlabel("r")
plt.ylabel("p_r")
plt.title("Ritratto di fase (mappa Poincaré)")
plt.show()
No description has been provided for this image
No description has been provided for this image
In [ ]: