Canonical Perturbation Theory

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

La Canonical Perturbation Theory studia sistemi hamiltoniani che differiscono debolmente da un sistema integrabile esatto. Si considera un Hamiltoniano della forma $$\large H(q,p) = H_0(q,p) + \varepsilon H_1(q,p), \qquad 0 < \varepsilon \ll 1, $$ dove $H_0$ è integrabile e $H_1$ rappresenta una perturbazione debole.

Il primo passo consiste nel passare alle coordinate azione--angolo $(I,\theta)$ associate a $H_0$, nelle quali $$\large H_0 = H_0(I), \qquad \dot{\theta} = \omega(I), \qquad \dot{I} = 0. $$ In queste variabili il moto non perturbato è una rotazione uniforme sul toro invariabile $I=\mathrm{cost}$.

La perturbazione introduce una dipendenza angolare: $$\large H(I,\theta) = H_0(I) + \varepsilon H_1(I,\theta), $$ che in generale rompe l’integrabilità esatta. L’obiettivo della teoria perturbativa canonica è costruire una trasformazione canonica quasi-identità $$\large (I,\theta) \mapsto (I',\theta') $$

tale da eliminare, ordine per ordine in $\varepsilon$, la dipendenza angolare non risonante dell’Hamiltoniano.

La trasformazione è generata da una funzione generatrice espansa in serie perturbativa: $$\large S(I',\theta) = I'\theta + \varepsilon S_1(I',\theta)\varepsilon^2 S_2(I',\theta) + \cdots $$

e conduce a un Hamiltoniano efficace $$\large H'(I') = H_0(I') + \varepsilon \langle H_1 \rangle\mathcal{O}(\varepsilon^2), $$ dove $\langle \cdot \rangle$ denota la media angolare.

In questo quadro, le azioni risultano invarianti fino a ordini elevati, mentre la frequenza del moto viene rinormalizzata: $$\large \langle \dot I \rangle = 0 + \mathcal{O}(\varepsilon^{n+1}). $$

come dire che $\textbf{I}=cost$ cioè $\textbf{I}$ rimane invatiato

$$\large \omega(I) \;\longrightarrow\; \omega(I) + \varepsilon\,\Delta\omega(I). $$

Questo fenomeno è alla base del metodo di Lindstedt--Poincaré, che elimina le divergenze secolari ridefinendo il tempo canonico.

Per sistemi a un grado di libertà, la dinamica rimane integrabile e le orbite nello spazio delle fasi sono curve di livello dell’Hamiltoniano. La perturbazione non induce caos, ma si manifesta come una modifica della relazione azione--frequenza. Nei sistemi con più gradi di libertà, invece, la teoria perturbativa rivela la persistenza parziale dei tori invarianti (KAM) e l’emergere di regioni risonanti e caotiche.

La Canonical Perturbation Theory fornisce quindi un ponte concettuale tra l’integrabilità esatta e la complessità della dinamica hamiltoniana reale, conservando la struttura geometrica fondamentale del flusso canonico.

Perturbazione - soluzione simbolica

import sympy as sp

def generating_function_expansion(q, P, t, epsilon):
    '''
    Creates an expansion of the generating function up to second
    order.
        :param q: Generalized coordinate.
        :param P: Generalized momentum.
        :param t: Time.
        :param epsilon: Perturbation parameter.
    :return: 
        Expanded generating function.
    '''
    F2_0 = q * P # Unperturbed generating function
    F2_1 = epsilon * (q ** 2 + P * t) # First-order correction example
    F2_2 = epsilon**2 * (sp.sin(q * sp.pi) * sp.cos(P * sp.pi)) # Second-order correction example
    return F2_0 + F2_1 + F2_2

q, P, t, epsilon = sp.symbols('q P t epsilon')
F2_expanded = generating_function_expansion(q, P, t, epsilon)

def hamilton_jacobi_expansion(q, P, t, epsilon):
    '''
    Solves the Hamilton-Jacobi equation using a power series
    expansion in epsilon.
        :param q: Generalized coordinate.
        :param P: Generalized momentum.
        :param t: Time.
        :param epsilon: Perturbation parameter.
    :return: 
        Action function series.
    '''
    S0	= P * q - 0.5 * P**2 * t # Example of unperturbed action
    S1	= epsilon * (P * sp.sin(q)) # First-order perturbation
    S2	= epsilon**2 * (P * sp.cos(q) * t) # Second-ordermperturbation
    return S0 + S1 + S2

S_expanded = hamilton_jacobi_expansion(q, P, t, epsilon) 

def lindstedt_poincare_method(omega_0, epsilon):
    '''
    Adjusts system's angular frequency using Lindstedt-Poincare.
        :param omega_0: Original angular frequency.
        :param epsilon: Perturbation parameter.
    :return: 
        Expanded angular frequency.
    '''
    omega_1 = epsilon * omega_0 / 10 # Example of first-order correction
    omega_2 = epsilon**2 * omega_0 / 100 # Example of second-order correction
    return omega_0 + omega_1 + omega_2

omega = lindstedt_poincare_method(1.0, epsilon)

def perturbed_harmonic_oscillator(m, omega_0, q, p, epsilon):
    '''
    Analyzes a perturbed harmonic oscillator.
        :param m: Mass of oscillator.
        :param omega_0: Natural frequency.
        :param q: Generalized coordinate.
        :param p: Generalized momentum.
        :param epsilon: Small perturbation parameter.
    :return: 
        Hamiltonian with perturbation.
    '''
    V = epsilon * sp.sin(q) # Example perturbative potential
    H0 = p**2 / (2*m)+0.5*m* omega_0**2 * q**2
    return H0 + V 
    
m, q, p = sp.symbols('m q p')

H_perturbed = perturbed_harmonic_oscillator(m, 1.0, q, p, epsilon)
print("Expanded Generating Function:")
display(F2_expanded)
print("Expanded Action Function from Hamilton-Jacobi:")
display(S_expanded)
print("Adjusted Frequency using Lindstedt-Poincare:")
display(omega)
print("Perturbed Harmonic Oscillator Hamiltonian:")
display(H_perturbed)

Expanded Generating Function:

$\displaystyle P q + \epsilon^{2} \sin{\left(\pi q \right)} \cos{\left(\pi P \right)} + \epsilon \left(P t + q^{2}\right)$

Expanded Action Function from Hamilton-Jacobi:

$\displaystyle - 0.5 P^{2} t + P \epsilon^{2} t \cos{\left(q \right)} + P \epsilon \sin{\left(q \right)} + P q$

Adjusted Frequency using Lindstedt-Poincare:

$\displaystyle 0.01 \epsilon^{2} + 0.1 \epsilon + 1.0$

Perturbed Harmonic Oscillator Hamiltonian:

$\displaystyle \epsilon \sin{\left(q \right)} + 0.5 m q^{2} + \frac{p^{2}}{2 m}$

Perturbazione - soluzione numerica - non derive secolari $\quad \rightarrow \quad$ nel piano delle fasi orbita chiusa

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

# -----------------------------
# PARAMETRI (tuning visivo)
# -----------------------------
epsilon = 1.65      # perturbazione ben visibile
Tmax = 1900.0        # tempo lungo -> precessione evidente
N = 60000           # alta risoluzione

# condizioni iniziali
q0 = 1.5
p0 = 0.0

# -----------------------------
# Equazioni di Hamilton
# -----------------------------
def hamiltonian_perturbed(t, y):
    q, p = y
    return [p, -q - epsilon * np.cos(q)]

def hamiltonian_unperturbed(t, y):
    q, p = y
    return [p, -q]

# -----------------------------
# Integrazione
# -----------------------------
t_eval = np.linspace(0, Tmax, N)

sol_p = solve_ivp(
    hamiltonian_perturbed,
    (0, Tmax), [q0, p0],
    t_eval=t_eval,
    rtol=1e-10, atol=1e-10
)

sol_u = solve_ivp(
    hamiltonian_unperturbed,
    (0, Tmax), [q0, p0],
    t_eval=t_eval,
    rtol=1e-10, atol=1e-10
)

q_p, p_p = sol_p.y
q_u, p_u = sol_u.y

# -----------------------------
# PLOT spazio delle fasi
# -----------------------------
plt.figure(figsize=(7,7))
for k in range(0, N, 4000):
    plt.plot(q_p[k:k+4000], p_p[k:k+4000], color='darkred', alpha=0.5)

plt.xlabel(r'$q$')
plt.ylabel(r'$p$')
plt.title('Precessione vista per cicli successivi')
plt.axis('equal')
plt.grid(alpha=0.3)
plt.show()

Perturbazioni - Mappa di Poincaré - circa un segmento

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

# -----------------------------
# Parametri
# -----------------------------
epsilon = 0.25
q0, p0 = 1.5, 0.0
Tmax = 800.0

# -----------------------------
# Equazioni di Hamilton
# -----------------------------
def hamiltonian(t, y):
    q, p = y
    return [p, -q - epsilon*np.cos(q)]

# -----------------------------
# Integrazione
# -----------------------------
t_eval = np.linspace(0, Tmax, 400000)

sol = solve_ivp(
    hamiltonian,
    (0, Tmax),
    [q0, p0],
    t_eval=t_eval,
    rtol=1e-10,
    atol=1e-10
)

q, p = sol.y

# -----------------------------
# Sezione di Poincaré: p = 0, dp/dt < 0
# -----------------------------
q_sec = []

for i in range(1, len(p)):
    if p[i-1] > 0 and p[i] <= 0:
        q_sec.append(q[i])

q_sec = np.array(q_sec)

# -----------------------------
# Plot
# -----------------------------
plt.figure(figsize=(6,4))
# non modifico la realtà perché devo guardare nel piano p=0
plt.scatter(q_sec, np.zeros_like(q_sec), s=10, color='darkred')

plt.xlabel(r'$q$')
plt.yticks([])
plt.title(r'Sezione di Poincaré: $p=0,\ \dot p<0$')
plt.grid(alpha=0.3)

plt.show()