KAM - Kolgomorov-Arnold-Moser Theorem¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
Consideriamo un sistema hamiltoniano integrabile¶
$$\large H_0(I), \quad I \in \mathbb{R}^n, $$con coordinate d'azione-angle $(I, \theta)$. Le traiettorie del sistema si muovono su tori invarianti con frequenze $$\large \omega(I) = \nabla H_0(I). $$ Se il sistema integrabile viene perturbato leggermente $$\large H(I, \theta) = H_0(I) + \epsilon H_1(I, \theta), \quad 0<\epsilon \ll 1, $$ il teorema KAM afferma che, sotto una condizione di non risonanza sulle frequenze $\omega(I)$ e alcune ipotesi di regolarità, una gran parte dei tori invarianti del sistema integrabile originale sopravvive, sebbene deformata dalla perturbazione. Questi tori sono detti $\textbf{tori KAM}$.
La condizione di non risonanza si scrive spesso come:¶
$$\large |\langle k, \omega(I) \rangle| \ge \frac{C}{|k|^\tau}, \quad \forall k \in \mathbb{Z}^n \setminus \{0\}, $$dove:
- $k = (k_1, \dots, k_n)$ è un vettore di numeri interi;
- $\omega(I) = (\omega_1(I), \dots, \omega_n(I))$ è il vettore delle frequenze del sistema integrabile;
- $\langle k, \omega(I) \rangle = k_1 \omega_1(I) + \dots + k_n \omega_n(I)$ è il \emph{prodotto scalare} tra $k$ e $\omega(I)$;
- $C>0$ e $\tau>n-1$ sono costanti che dipendono dal sistema.
$\textbf{Significato del primo membro:}$ $$\large |\langle k, \omega(I) \rangle| $$ rappresenta la distanza (in valore assoluto) della combinazione lineare delle frequenze pesata dagli interi $k$ dallo zero.
In altre parole, misura quanto la combinazione lineare $$\large k_1 \omega_1 + \dots + k_n \omega_n $$ si avvicina a una risonanza $(=0)$.
- Se fosse esattamente zero ($\langle k, \omega(I) \rangle = 0$) per qualche $k \neq 0$, le frequenze sarebbero \emph{risonanti}, e il moto sul toro potrebbe chiudersi su se stesso dopo un certo numero di cicli. In questo caso, anche piccole perturbazioni possono distruggere il toro.
- La condizione di non risonanza assicura che nessuna combinazione lineare intera non banale delle frequenze sia troppo vicina a zero, come quantificato dal lato destro $\frac{C}{|k|^\tau}$.
In pratica, il primo membro misura quanto le frequenze \emph{evitano le risonanze}, condizione fondamentale per la persistenza dei tori invarianti sotto piccole perturbazioni.
Esercizio del libro¶
Confronta le vecchie Azioni - poco significativo¶
In [79]:
import numpy as np
def hamiltonian(I, theta, H0, H1, epsilon):
'''
Calculate the perturbed Hamiltonian.
:param I: Action variables.
:param theta: Angle variables.
:param HO: Unperturbed Hamiltonian function.
:param Hl: Perturbation Hamiltonian function.
:param epsilon: Small perturbation parameter.
:return:
Total Hamiltonian.
'''
return H0(I, theta) + epsilon * H1(I, theta)
def validate_non_resonance(omega, tau, C):
'''
Validate the non-resonance condition for frequencies.
:param omega: Frequency vector.
:param tau: Resonance condition parameter.
:param C: Positive constant for non-resonance.
:return:
True if non-resonance condition is met, False otherwise.
'''
omega1 = np.array(omega)
for k in np.ndindex(omega1.shape):
if np.linalg.norm(k) != 0: # Ensure k is not the zero vector
if abs(sum(k * omega1)) < C / (np.linalg.norm(k) ** tau):
return False
return True
def iterative_kam_validation(I_init, H0, epsilon, tau, C, iterations=3000):
'''
Iteratively validate invariant tori with KAM conditions.
:param I_init: Initial action variables.
:param HO: Unperturbed Hamiltonian function.
:param epsilon: Small perturbation parameter.
:param tau: Resonance condition parameter.
:param C: Positive constant for non-resonance.
:param iterations: Number of iterations for validation.
:return:
Final state of action variables and validation status.
'''
I = I_init
for _ in range(iterations):
omega = np.gradient(H0(I))
if not validate_non_resonance(omega, tau, C):
return I, False
I += -epsilon * np.gradient(epsilon * np.random.rand(*I.shape)) # Example update
return I, True
# Example functions for unperturbed and perturbed Hamiltonians
def H0_example(I):
return np.sum(I**2 / 2)
def H1_example(I, theta):
return np.sum(np.sin(theta) * I)
# Example usage
I_init = np.array([1.0, 2.0])
theta = 1
epsilon = 0.01
tau = 2.5
C = 1.0
final_I, is_valid = iterative_kam_validation(I_init, H0_example, epsilon, tau, C)
print("Final action variables:", final_I)
print("KAM validation status:", is_valid)
Final action variables: [0.99999408 1.99999408] KAM validation status: True
Toy Hamiltoniano Giove–Saturno - con due metodi + controllo finale¶
In [101]:
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.integrate import solve_ivp
# -----------------------------
# 1. Hamiltoniano Giove-Saturno toy
# -----------------------------
def H0(I):
# parte integrabile
return 0.5 * np.sum(I**2)
def H1(I, theta):
# perturbazione tipo near 5:2
return np.cos(5*theta[1] - 2*theta[0])
def equations(t, y, epsilon):
I = y[:2]
theta = y[2:]
# Frequenze integrabili
omega = I
# Derivate
dtheta = omega + epsilon * np.array([
2*np.sin(5*theta[1] - 2*theta[0]),
-5*np.sin(5*theta[1] - 2*theta[0])
])
dI = epsilon * np.array([
2*np.sin(5*theta[1] - 2*theta[0]),
-5*np.sin(5*theta[1] - 2*theta[0])
])
return np.concatenate([dI, dtheta])
# -----------------------------
# 2. Integrazione simplettica (Stormer-Verlet)
# -----------------------------
def symplectic_integrator(I0, theta0, epsilon, dt, steps):
I = I0.copy()
theta = theta0.copy()
I_hist = []
theta_hist = []
for _ in range(steps):
# kick
I_half = I + 0.5 * dt * epsilon * np.array([
2*np.sin(5*theta[1] - 2*theta[0]),
-5*np.sin(5*theta[1] - 2*theta[0])
])
# drift
theta = theta + dt * I_half
# kick
I = I_half + 0.5 * dt * epsilon * np.array([
2*np.sin(5*theta[1] - 2*theta[0]),
-5*np.sin(5*theta[1] - 2*theta[0])
])
I_hist.append(I.copy())
theta_hist.append(theta.copy())
return np.array(I_hist), np.array(theta_hist)
# -----------------------------
# 3. Parametri realistici Giove-Saturno
# -----------------------------
I0 = np.array([5.2, 9.5])
theta0 = np.array([0.1, 0.2])
epsilon = 1e-3
dt = 0.01
steps = 200000 # tempo lungo ma non folle
# -----------------------------
# 4. Simulazione
# -----------------------------
#I_hist = symplectic_integrator(I0, theta0, epsilon, dt, steps)
I_hist, theta_hist = symplectic_integrator(
I0, theta0, epsilon, dt, steps
)
# -----------------------------
# 5. Diagnostica KAM-like
# -----------------------------
# Oscillazione massima
delta_I = np.max(np.linalg.norm(I_hist - I0, axis=1))
print("Max |I(t) - I(0)| =", delta_I)
# Azione mediata
I_mean = np.mean(I_hist, axis=0)
print("Time-averaged action:", I_mean)
# -----------------------------
# 6. Frequency Map Analysis minimale
# -----------------------------
def dominant_frequency(signal, dt):
fft = np.abs(rfft(signal - np.mean(signal)))
freqs = rfftfreq(len(signal), dt)
return freqs[np.argmax(fft[1:]) + 1]
def fundamental_frequency_from_slope(theta, dt):
t = np.arange(len(theta)) * dt
coeff = np.polyfit(t, theta, 1)
return coeff[0] # pendenza = frequenza
from numpy.fft import rfft, rfftfreq
from numpy.fft import fft, fftfreq
def fundamental_frequency_fft(theta, dt):
z = np.exp(1j * theta)
z = z - np.mean(z)
F = np.abs(fft(z))
freqs = fftfreq(len(z), dt)
# solo frequenze positive
mask = freqs > 0
return freqs[mask][np.argmax(F[mask])]
# Metodo A — Frequenze come pendenza media (molto robusto)
omega1 = fundamental_frequency_from_slope(theta_hist[:,0], dt)
omega2 = fundamental_frequency_from_slope(theta_hist[:,1], dt)
print("Fundamental frequencies (slope):", omega1, omega2)
print("Ratio omega2 / omega1:", omega2 / omega1)
print('')
# Metodo B — Frequency Map (FFT su e^(i*theta) )
omega1_fft = fundamental_frequency_fft(theta_hist[:,0], dt)
omega2_fft = fundamental_frequency_fft(theta_hist[:,1], dt)
print("Fundamental frequencies (FFT):", omega1_fft, omega2_fft)
print("Ratio omega2 / omega1:", omega2_fft / omega1_fft)
print('')
# Ripeto l’estrazione delle frequenze su due finestre temporali e confronto
mid = len(theta_hist) // 2
omega1_a = fundamental_frequency_from_slope(theta_hist[:mid,0], dt)
omega1_b = fundamental_frequency_from_slope(theta_hist[mid:,0], dt)
print("Frequency drift:", abs(omega1_a - omega1_b))
Max |I(t) - I(0)| = 0.0002434546018488399 Time-averaged action: [5.20003713 9.49990718] Fundamental frequencies (slope): 5.200037127652402 9.499907180868544 Ratio omega2 / omega1: 1.8268921832020366 Fundamental frequencies (FFT): 0.8275 1.512 Ratio omega2 / omega1: 1.827190332326284 Frequency drift: 1.4033219031261979e-13
In [ ]: