Chaos in Hamiltonian Systems

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

Il modello di Hénon--Heiles

Il modello di Hénon--Heiles fu introdotto nel 1964 da Michel Hénon e Carl Heiles con l'obiettivo di studiare la stabilità delle orbite stellari in potenziali galattici debolmente non integrabili.

L'idea fondamentale è quella di considerare un sistema hamiltoniano conservativo a due gradi di libertà che rappresenti il moto di una stella nel piano di una galassia non perfettamente simmetrica. Il sistema di partenza è un oscillatore armonico bidimensionale, che è completamente integrabile, al quale viene aggiunta una perturbazione non lineare minimale che rompe l'integrabilità ma conserva l'energia.

Nel modello, le coordinate canoniche $(q_1, q_2)$ rappresentano la posizione della stella nel piano galattico, mentre i momenti coniugati $(p_1, p_2)$ rappresentano le quantità di moto associate. Il potenziale risultante descrive una galassia quasi ellittica, con una debole triaxialità o asimmetria.

Un aspetto cruciale del modello è che, pur essendo completamente deterministico e privo di dissipazione o forzanti esterne, esso può esibire caos deterministico. Il comportamento dinamico dipende in modo essenziale dall'energia totale del sistema.

Per basse energie, il moto è regolare e le traiettorie nello spazio delle fasi sono confinate su tori invarianti, in accordo con il teorema KAM. All'aumentare dell'energia, i tori iniziano a rompersi, compaiono risonanze e regioni caotiche, e la dinamica diventa progressivamente più complessa. Per energie sufficientemente elevate, alcune traiettorie possono addirittura sfuggire dal potenziale.

Il modello di Hénon--Heiles ha avuto un ruolo storico fondamentale perché ha mostrato per la prima volta in modo chiaro che il caos può emergere in sistemi hamiltoniani conservativi. Ancora oggi viene utilizzato come sistema prototipo per lo studio delle sezioni di Poincaré, degli esponenti di Lyapunov e della transizione ordine--caos nei sistemi dinamici non lineari.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def hamiltonian(q, p, epsilon):
    q1, q2 = q
    p1, p2 = p

    kinetic = 0.5 * (p1**2 + p2**2)
    potential = 0.5 * (q1**2 + q2**2) \
                + epsilon * (q1**2 * q2 - q2**3 / 3)

    return kinetic + potential

def hamiltonian_equations(t, y, epsilon):
    q1, q2, p1, p2 = y

    dq1 = p1
    dq2 = p2

    dp1 = -q1 - 2 * epsilon * q1 * q2
    dp2 = -q2 - epsilon * (q1**2 - q2**2)

    return [dq1, dq2, dp1, dp2]

def lyapunov_exponent(y0, t_span, perturbation, epsilon, dt=0.01):
    y = np.array(y0)
    y_pert = y + perturbation

    t0, tf = t_span
    times = np.arange(t0, tf, dt)

    delta0 = np.linalg.norm(perturbation)
    lyap_sum = 0.0

    for _ in times:
        sol1 = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y,
            args=(epsilon,),
            rtol=1e-9
        )

        sol2 = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y_pert,
            args=(epsilon,),
            rtol=1e-9
        )

        y = sol1.y[:, -1]
        y_pert = sol2.y[:, -1]

        delta = y_pert - y
        dist = np.linalg.norm(delta)

        lyap_sum += np.log(dist / delta0)

        delta = delta / dist * delta0
        y_pert = y + delta

    return lyap_sum / (tf - t0)

def poincare_section(y0, p_section, t_span, epsilon, dt=0.01):
    t0, tf = t_span
    times = np.arange(t0, tf, dt)

    y = np.array(y0)
    points = []

    q1_prev = y[0]

    for _ in times:
        sol = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y,
            args=(epsilon,),
            rtol=1e-9
        )

        y = sol.y[:, -1]
        q1 = y[0]

        if q1_prev < p_section and q1 >= p_section and y[2] > 0:
            points.append([y[1], y[3]])

        q1_prev = q1

    return np.array(points)

epsilon = 1.0
t_span = (0, 3000)

y0 = [0.0, 0.10, 0.30, 0.0] # tosus liscio - Energia circa 0.055
y0 = [0.0, 0.20, 0.45, 0.0] # isole KAM - Energia circa 0.095
#y0 = [0.0, 0.25, 0.50, 0.0] # caos con isole immerse - Energia circa 0.125


points = poincare_section(y0, 0.0, t_span, epsilon)

plt.figure(figsize=(5, 5))
plt.scatter(points[:, 0], points[:, 1], s=1)
plt.xlabel("q2")
plt.ylabel("p2")
plt.title("Sezione di Poincaré (q1 = 0)")
plt.show()

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def hamiltonian(q, p, epsilon):
    q1, q2 = q
    p1, p2 = p

    kinetic = 0.5 * (p1**2 + p2**2)
    potential = 0.5 * (q1**2 + q2**2) \
                + epsilon * (q1**2 * q2 - q2**3 / 3)

    return kinetic + potential

def hamiltonian_equations(t, y, epsilon):
    q1, q2, p1, p2 = y

    dq1 = p1
    dq2 = p2

    dp1 = -q1 - 2 * epsilon * q1 * q2
    dp2 = -q2 - epsilon * (q1**2 - q2**2)

    return [dq1, dq2, dp1, dp2]

def lyapunov_exponent(y0, t_span, perturbation, epsilon, dt=0.01):
    y = np.array(y0)
    y_pert = y + perturbation

    t0, tf = t_span
    times = np.arange(t0, tf, dt)

    delta0 = np.linalg.norm(perturbation)
    lyap_sum = 0.0

    for _ in times:
        sol1 = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y,
            args=(epsilon,),
            rtol=1e-9
        )

        sol2 = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y_pert,
            args=(epsilon,),
            rtol=1e-9
        )

        y = sol1.y[:, -1]
        y_pert = sol2.y[:, -1]

        delta = y_pert - y
        dist = np.linalg.norm(delta)

        lyap_sum += np.log(dist / delta0)

        delta = delta / dist * delta0
        y_pert = y + delta

    return lyap_sum / (tf - t0)

def poincare_section(y0, p_section, t_span, epsilon, dt=0.01):
    t0, tf = t_span
    times = np.arange(t0, tf, dt)

    y = np.array(y0)
    points = []

    q1_prev = y[0]

    for _ in times:
        sol = solve_ivp(
            hamiltonian_equations,
            [0, dt],
            y,
            args=(epsilon,),
            rtol=1e-9
        )

        y = sol.y[:, -1]
        q1 = y[0]

        if q1_prev < p_section and q1 >= p_section and y[2] > 0:
            points.append([y[1], y[3]])

        q1_prev = q1

    return np.array(points)

epsilon = 1.0
t_span = (0, 3000)

y0 = [0.0, 0.10, 0.30, 0.0] # tosus liscio - Energia circa 0.055
y0 = [0.0, 0.20, 0.45, 0.0] # isole KAM - Energia circa 0.095
y0 = [0.0, 0.25, 0.50, 0.0] # caos con isole immerse - Energia circa 0.125


points = poincare_section(y0, 0.0, t_span, epsilon)

plt.figure(figsize=(5, 5))
plt.scatter(points[:, 0], points[:, 1], s=1)
plt.xlabel("q2")
plt.ylabel("p2")
plt.title("Sezione di Poincaré (q1 = 0)")
plt.show()

Sì: l’attrito mareale è onnipresente nel Sistema Solare reale e contribuisce alla stabilità.

• Dissipa energia
• Smorza eccentricità e inclinazioni
• Porta verso stati “ordinati” (risonanze, sincronizzazioni)

Esempi concreti:

• sincronizzazione Luna–Terra
• circolarizzazione delle orbite dei satelliti
• cattura in risonanza (es. Io–Europa–Ganimede)
• Quanto detto non è caos hamiltoniano puro.

️Perché allora si parla tanto di KAM?

Il teorema KAM non pretende di descrivere il Sistema Solare reale in toto. Serve a rispondere a una domanda più limitata:

• se ignorassi completamente la dissipazione, il sistema sarebbe catastroficamente instabile?
• La risposta (sorprendente) è: no, una grande frazione delle orbite resta quasi-periodica.

KAM è quindi:

• un baseline teorico,
• un limite ideale,
• non un modello completo della realtà.

Attrito mareale: amico della stabilità, nemico di KAM

Qui la tua intuizione è perfetta.

• Il teorema KAM richiede:
• sistema conservativo
• hamiltoniano
• perturbazioni piccole

L’attrito mareale:

• rompe la struttura simplettica
• distrugge i tori invarianti
• introduce attrattori

Quindi:

• KAM non è applicabile in senso stretto a un sistema con attrito.
• Ma questo non è un difetto di KAM: è un cambio di paradigma dinamico.

Paradosso solo apparente

Sembra paradossale:

• KAM dice che il sistema è stabile senza dissipazione

tu osservi che con dissipazione è ancora più stabile

La risoluzione è:

• KAM parla di stabilità neutra (tori)
• la dissipazione introduce stabilità asintotica

Sono due nozioni diverse di stabilità.

Come si modellizza davvero il Sistema Solare. In pratica si fa questo:

Passo 1

— Sistema hamiltoniano quasi-integrabile

• problema dei N corpi
• analisi KAM / risonanze
• si capisce dove può nascere il caos

Passo 2 — Correzioni dissipative secolari

• maree
• relatività generale
• effetti non conservativi
• Questi termini agiscono su tempi lunghissimi
• non creano il caos
• ma guidano lentamente l’evoluzione

Il caos determina dove puoi andare, la dissipazione determina dove finisci.

Visione moderna (molto importante)

Il Sistema Solare è dinamicamente “marginale”: quasi integrabile, leggermente caotico, debolmente dissipativo.

In questo regime:

• KAM è rilevante per la struttura dello spazio delle fasi
• la dissipazione governa l’evoluzione secolare

Un’analogia utile: un tavolo quasi piano

• KAM descrive le valli e le creste del paesaggio
• l’attrito fa sì che una pallina non rimbalzi all’infinito, ma scivoli lentamente verso una valle

Senza KAM:

• il paesaggio sarebbe caotico ovunque, la pallina cadrebbe subito

È sensato modellizzare il Sistema Solare con KAM se c’è attrito mareale?

• Sì, come struttura portante dello spazio delle fasi
• No, come descrizione completa dell’evoluzione temporale
• l’attrito è fondamentale per la stabilità reale,
• KAM è fondamentale per capire perché il sistema non esplode caoticamente anche senza attrito.

In una frase conclusiva

• KAM spiega perché il Sistema Solare può esistere; la dissipazione spiega perché persiste.