Hamiltonian Chaos and Lyapunov Exponents¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
Caos hamiltoniano ed esponenti di Lyapunov nel doppio pendolo¶
Il doppio pendolo costituisce uno dei sistemi paradigmatici per lo studio del $\textit{caos hamiltoniano}$. Esso è un sistema conservativo a due gradi di libertà, non integrabile in forma chiusa, che presenta una ricca struttura dello spazio delle fasi caratterizzata dalla coesistenza di regioni regolari e caotiche.
Modello dinamico
Il sistema può essere descritto tramite le coordinate generalizzate $(\theta_1, \theta_2)$ e le rispettive velocità angolari $(\dot{\theta}_1, \dot{\theta}_2)$. In assenza di dissipazione, la dinamica deriva da un Hamiltoniano del tipo $$\large H = T(\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2) + V(\theta_1,\theta_2), $$ che garantisce la conservazione dell’energia totale.
Per piccole oscillazioni attorno all’equilibrio stabile $(\theta_1,\theta_2) = (0,0)$ il sistema è approssimabile da un modello quasi-lineare e risulta localmente integrabile.
Stabilità e regime quasi-integrabile
Nel limite di basse energie il moto è dominato dai modi normali del sistema (sincrono e antisincrono). In questa regione lo spazio delle fasi è foliato da tori invarianti, in accordo con il teorema di Kolmogorov–Arnold–Moser (KAM). Le traiettorie risultano stabili e la separazione di due orbite inizialmente vicine cresce al più in modo polinomiale nel tempo.
Il massimo esponente di Lyapunov $$\large \lambda_1 = \lim_{t \to \infty} \frac{1}{t} \ln \frac{\|\delta \mathbf{x}(t)\|}{\|\delta \mathbf{x}(0)\|} $$ risulta nullo entro l’errore numerico: $$\large \lambda_1 \approx 0. $$
Transizione al caos
All’aumentare dell’energia (ovvero dell’ampiezza iniziale degli angoli), le risonanze tra i modi normali diventano rilevanti. Per valori tipici degli angoli iniziali $$\large \theta \simeq 1.2\text{–}1.3 \ \text{rad}, $$ si osserva la distruzione progressiva dei tori KAM e la nascita di regioni iperboliche nello spazio delle fasi.
In questa regione la dinamica diventa altamente sensibile alle condizioni iniziali e il massimo esponente di Lyapunov assume valori positivi: $$\large \lambda_1 > 0, $$ segnalando la presenza di caos deterministico.
Spazio delle fasi misto
Anche nel regime caotico, il sistema non è globalmente ergodico. Lo spazio delle fasi presenta una struttura mista, con isole di stabilità immerse in un mare caotico. Il valore di $\lambda_1$ dipende quindi dalla traiettoria considerata e dall’energia iniziale.
Aspetti numerici
Il calcolo corretto degli esponenti di Lyapunov in sistemi hamiltoniani richiede:
- equazioni del moto corrette e conservative;
- integratori numerici stabili (preferibilmente simplettici);
- rinormalizzazione periodica della perturbazione;
- tempi di integrazione sufficientemente lunghi.
Un valore positivo e robusto di $\lambda_1$, indipendente dal passo temporale e dall’ampiezza iniziale della perturbazione, costituisce una chiara firma del caos hamiltoniano.
import numpy as np
# -----------------------------
# Parametri fisici
# -----------------------------
g = 9.81
L = 1.0
m = 1.0
# -----------------------------
# Equazioni del doppio pendolo
# -----------------------------
def double_pendulum(t, y):
theta1, omega1, theta2, omega2 = y
delta = theta2 - theta1
den = 2 - np.cos(delta)**2
dtheta1 = omega1
dtheta2 = omega2
domega1 = (
-2*np.sin(theta1)
+ np.sin(theta2)*np.cos(delta)
- omega2**2*np.sin(delta)
- omega1**2*np.sin(delta)*np.cos(delta)
) / den
domega2 = (
2*np.sin(theta1)*np.cos(delta)
- 2*np.sin(theta2)
+ 2*omega1**2*np.sin(delta)
+ omega2**2*np.sin(delta)*np.cos(delta)
) / den
return np.array([dtheta1, domega1, dtheta2, domega2])
# -----------------------------
# Integratore RK4
# -----------------------------
def rk4_step(f, t, y, dt):
k1 = f(t, y)
k2 = f(t + dt/2, y + dt*k1/2)
k3 = f(t + dt/2, y + dt*k2/2)
k4 = f(t + dt, y + dt*k3)
return y + dt * (k1 + 2*k2 + 2*k3 + k4) / 6
# -----------------------------
# Calcolo esponente di Lyapunov
# -----------------------------
def lyapunov_exponent(
y0,
delta0=1e-8,
dt=0.01,
n_steps=200000,
renorm_steps=10
):
y = np.array(y0, dtype=float)
# y_pert = y + delta0 * np.random.randn(len(y0))
y_pert = y.copy()
y_pert[0] += delta0
t = 0.0
sum_log = 0.0
count = 0
for i in range(n_steps):
y = rk4_step(double_pendulum, t, y, dt)
y_pert = rk4_step(double_pendulum, t, y_pert, dt)
t += dt
if i % renorm_steps == 0:
delta_vec = y_pert - y
dist = np.linalg.norm(delta_vec)
if dist == 0:
continue
sum_log += np.log(dist / delta0)
count += 1
# rinormalizzazione
delta_vec = delta_vec / dist * delta0
y_pert = y + delta_vec
return sum_log / (count * dt * renorm_steps)
# -----------------------------
# Condizioni iniziali caotiche
# -----------------------------
icount = 5
if icount == 1:
# condizioni iniziali caotiche
# Massimo esponente di Lyapunov λ₁ ≈ 0.1674
y0 = [
np.pi / 2, # theta1
0.0, # omega1
np.pi / 2, # theta2
0.0 # omega2
]
if icount == 2:
# Punto fisso stabile (oscillazioni piccole)
# Massimo esponente di Lyapunov λ₁ ≈ 0.0008
y0 = [
0.05, # theta1 (rad)
0.0, # omega1
0.05, # theta2
0.0 # omega2
]
if icount == 3:
# Modo normale in fase (stabile)
# Massimo esponente di Lyapunov λ₁ ≈ 0.0024
# Qui stai eccitando il modo normale simmetrico
y0 = [
0.3, # theta1
0.0,
0.3, # theta2
0.0
]
if icount == 4:
# Modo normale fuori fase (ancora stabile a bassa energia)
# Questo è il modo antisimmetrico
# Massimo esponente di Lyapunov λ₁ ≈ 0.0037
y0 = [
0.3, # theta1
0.0,
-0.3, # theta2
0.0
]
if icount == 5:
# Modo normale fuori fase (non stabile)
# Questo è il modo antisimmetrico
# Massimo esponente di Lyapunov λ₁ ≈ 0.1722
# per angoli = 1.2, -1,2 è stabile
# a 1.3, -1.3 scatta improvvisamente a instabilità
y0 = [
1.3, # theta1
0.0,
-1.3, # theta2
0.0
]
lambda1 = lyapunov_exponent(y0)
print(f"Massimo esponente di Lyapunov λ₁ ≈ {lambda1:.4f}")
# l'esecuzione richiede circa 30 secondi
Massimo esponente di Lyapunov λ₁ ≈ 0.1722