Phase Space and Distribution Functions¶

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

Osillatore armonico¶

Il grafico p, q mostra una ellisse stabile --> H = costante = conservazione del volume nella evoluzione temporale¶

In [1]:
# Phase space & distribution functions example
import numpy as np
import matplotlib.pyplot as plt

# Hamiltonian
def H(q, p):
    return 0.5 * (p**2 + q**2)

# Hamilton's equations
def hamilton_equations(q, p, t):
    dqdt = p
    dpdt = -q
    return dqdt, dpdt

# Equations of motion (integrator: symplectic Euler)
def equations_of_motion(q0, p0, t):
    q = np.zeros_like(t)
    p = np.zeros_like(t)
    q[0], p[0] = q0, p0
    dt = t[1] - t[0]
    for i in range(len(t)-1):
        p[i+1] = p[i] - q[i] * dt
        q[i+1] = q[i] + p[i+1] * dt
    return q, p

# Liouville equation along characteristics
def liouville_equation_characteristics(rho0, q0, p0, t):
    q, p = equations_of_motion(q0, p0, t)
    rho = rho0 * np.ones_like(t)  # conserved along trajectories
    return q, p, rho

# Equilibrium distribution (canonical ensemble)
def equilibrium_distribution(q, p, beta):
    return np.exp(-beta * H(q, p))

# Time grid
t = np.linspace(0, 20, 2000)

# Initial conditions
q0, p0 = 1.0, 0.0
rho0 = 1.0
beta = 1.0

# Evolution
q, p = equations_of_motion(q0, p0, t)

# Plot: phase space trajectory
plt.figure()
plt.plot(q, p)
plt.xlabel("q")
plt.ylabel("p")
plt.title("Phase space trajectory (harmonic oscillator)")
plt.show()
No description has been provided for this image

Evoluzione temporale di una distribuzione di fase $\rho(q,p,t)$ per oscillatore armonico¶

Consideriamo un sistema hamiltoniano unidimensionale con Hamiltoniana $$\large H(q,p) = \frac{1}{2}(p^2 + q^2), $$ cioè l’oscillatore armonico classico.

Lo spazio delle fasi è bidimensionale, con coordinate $(q,p)$, e la dinamica è governata dalle equazioni di Hamilton $$\large \dot q = \frac{\partial H}{\partial p} = p, \qquad \dot p = -\frac{\partial H}{\partial q} = -q. $$

Distribuzione iniziale¶

La distribuzione iniziale $\rho(q,p,0)$ viene campionata come una gaussiana centrata nell’origine: $$\large \rho(q,p,0) \propto \exp\!\left(-\frac{q^2 + p^2}{2\sigma^2}\right). $$ Operativamente, ciò equivale a generare una nuvola di punti nel phase space, ciascuno dei quali rappresenta una caratteristica della dinamica.

Equazione di Liouville¶

L’evoluzione temporale della densità di fase è descritta dall’equazione di Liouville: $$\large \frac{\partial \rho}{\partial t}+ \dot q \frac{\partial \rho}{\partial q} + \dot p \frac{\partial \rho}{\partial p}=0$$

Lungo le traiettorie hamiltoniane $(q(t),p(t))$, questa equazione implica: $$\large \frac{d}{dt}\rho(q(t),p(t),t) = 0, $$ ossia la densità è conservata lungo le caratteristiche.

Metodo delle caratteristiche}¶

Nel metodo delle caratteristiche:

  • ogni punto iniziale $(q_0,p_0)$ evolve secondo le equazioni di Hamilton;
  • il valore di $\rho$ associato a ciascun punto resta costante;
  • la distribuzione evolve per trasporto nello spazio delle fasi.

Non si risolve quindi una PDE su una griglia, ma si segue l’evoluzione di un insieme di punti campionati.

Interpretazione geometrica}¶

Per l’oscillatore armonico, la soluzione esatta può essere scritta come una rotazione nel phase space: $$\large \begin{pmatrix} q(t) \\ p(t) \end{pmatrix} = \begin{pmatrix} \cos t & \sin t \\ -\sin t & \cos t \end{pmatrix} \begin{pmatrix} q_0 \\ p_0 \end{pmatrix}. $$

Ne consegue che:

  • le orbite sono chiuse (curve di energia costante);
  • la distribuzione $\rho(q,p,t)$ ruota rigidamente;
  • il volume di fase è conservato.

La conservazione del volume è espressione diretta del teorema di Liouville: $$\large \frac{d}{dt} \int_{\Omega(t)} dq\,dp = 0. $$

Assenza di rilassamento}¶

In questo sistema:

  • la dinamica è integrabile;
  • non vi è caos;
  • non vi è mixing nel phase space.

Pertanto, una distribuzione iniziale arbitraria non converge spontaneamente alla distribuzione canonica $$\large \rho_{\mathrm{eq}}(q,p) \propto e^{-\beta H(q,p)}. $$

La distribuzione resta semplicemente in rotazione, senza perdita di informazione fine-grained.

Osservazione fondamentale}¶

Questo esempio mostra chiaramente che:

  • l’equazione di Liouville non implica termalizzazione;
  • la conservazione della densità di fase non garantisce l’equilibrio statistico;
  • l’equilibrio richiede meccanismi aggiuntivi (caos, mixing, collisioni).

L’oscillatore armonico rappresenta quindi un caso didattico ideale per separare chiaramente dinamica microscopica ed equilibrio statistico.

In [3]:
# Evolution of a phase-space distribution rho(q,p,t)
import numpy as np
import matplotlib.pyplot as plt

# Hamiltonian and equations of motion (harmonic oscillator)
def equations_of_motion(q, p, dt):
    # symplectic Euler step
    p_new = p - q * dt
    q_new = q + p_new * dt
    return q_new, p_new

# Initial phase-space distribution: Gaussian blob
np.random.seed(0)
N = 5000
q = np.random.normal(0.0, 0.5, N)
p = np.random.normal(0.0, 0.5, N)

# Time parameters
dt = 0.05
steps = 200

# Store snapshots
snapshots = []
for n in [0, 50, 100, 150, 200]:
    qn, pn = q.copy(), p.copy()
    for _ in range(n):
        qn, pn = equations_of_motion(qn, pn, dt)
    snapshots.append((qn, pn))

# Plot snapshots
plt.figure()
for i, (qs, ps) in enumerate(snapshots):
    plt.scatter(qs, ps, s=1, alpha=0.5, label=f"t={i*50*dt:.2f}")
plt.xlabel("q")
plt.ylabel("p")
plt.title("Evolution of phase-space distribution ρ(q,p,t)")
plt.legend(markerscale=5)
plt.show()
No description has been provided for this image

Oscillatore anarmonico e distribuzione di fase¶

Consideriamo il sistema hamiltoniano unidimensionale con Hamiltoniana $$\large H(q,p) = \frac{p^2}{2} + \frac{q^4}{4}, $$ che descrive un oscillatore anarmonico con potenziale quartico.

Lo spazio delle fasi è bidimensionale e la dinamica è governata dalle equazioni di Hamilton: $$\large \dot q = \frac{\partial H}{\partial p} = p, \qquad \dot p = -\frac{\partial H}{\partial q} = -q^3. $$

Equazione di Liouville¶

La densità di fase $\rho(q,p,t)$ soddisfa l’equazione di Liouville: $$\large \frac{\partial \rho}{\partial t} p \frac{\partial \rho}{\partial q} q^3 \frac{\partial \rho}{\partial p}= 0. $$

Lungo le traiettorie hamiltoniane $(q(t),p(t))$ vale: $$\large \frac{d}{dt}\rho(q(t),p(t),t) = 0, $$ per cui la densità è conservata lungo le caratteristiche.

Distribuzione iniziale¶

Si considera una distribuzione iniziale gaussiana: $$\large \rho(q,p,0) \propto \exp\!\left(-\frac{q^2 + p^2}{2\sigma^2}\right), $$ campionata tramite una nuvola di punti nello spazio delle fasi.

Evoluzione nel phase space¶

A differenza dell’oscillatore armonico, la frequenza del moto dipende dall’energia. Ne consegue che:

  • le orbite di fase non sono circolari;
  • punti con energie diverse ruotano con velocità angolare diversa;
  • la distribuzione $\rho(q,p,t)$ subisce uno \emph{shearing} progressivo.

La densità fine-grained resta costante, ma la distribuzione sviluppa strutture filamentose sempre più sottili.

Osservazione fisica¶

Questo esempio mostra che:

  • la conservazione del volume di fase resta valida (teorema di Liouville);
  • l’anarmonicità introduce mixing cinematico;
  • il rilassamento osservabile è di tipo \emph{coarse-grained}, non microscopico.

L’oscillatore anarmonico rappresenta quindi un ponte concettuale tra sistemi integrabili e dinamica caotica.

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

# -----------------------------
# Hamilton equations (anharmonic oscillator)
# H = p^2/2 + q^4/4
# dq/dt = p
# dp/dt = -q^3
# -----------------------------
def equations_of_motion(q, p, dt):
    # simplectic Euler
    p_new = p - q**3 * dt
    q_new = q + p_new * dt
    return q_new, p_new

# -----------------------------
# Initial phase-space distribution
# -----------------------------
np.random.seed(0)
N = 5000
q0 = np.random.normal(0.0, 0.6, N)
p0 = np.random.normal(0.0, 0.6, N)

# -----------------------------
# Time parameters
# -----------------------------
dt = 0.02
times = [0.0, 2.0, 4.0, 6.0, 8.0]
steps = [int(t / dt) for t in times]

# -----------------------------
# Compute snapshots of ρ(q,p,t)
# -----------------------------
snapshots = []
for n in steps:
    q, p = q0.copy(), p0.copy()
    for _ in range(n):
        q, p = equations_of_motion(q, p, dt)
    snapshots.append((q, p))

# -----------------------------
# Plot phase-space evolution
# -----------------------------
plt.figure()
for (q, p), t in zip(snapshots, times):
    plt.scatter(q, p, s=1, alpha=0.4, label=f"t = {t}")
plt.xlabel("q")
plt.ylabel("p")
plt.title("Anharmonic oscillator: evolution of ρ(q,p,t)")
plt.legend(markerscale=5)
plt.show()
No description has been provided for this image
In [ ]:
In [7]:
import numpy as np 

def lagrangian(T, V):
    """
    Calculate the Lagrangian of a system.
        :param T: Kinetic energy.
        :param V: Potential energy.
    :return: 
        Lagrangian.
    """
    return T - V

def rayleigh_dissipation(velocities, coefficients):
    """
    Calculate the Rayleigh dissipation function.
    :param velocities: Generalized velocities.
    :param coefficients: Damping coefficients.
    :retum: 
        Rayleigh dissipation.
    """
    return 0.5 * np.sum(coefficients * velocities**2)

def harmonic_oscillator(m, k, gamma, x, v, dt):
    """
    Simulates a damped harmonic oscillator using the velocity Verlet method.
        :param m: Mass of the oscillator.
        :param k: Spring constant.
        :param gamma: Damping coefficient.
        :param x: Initial position.
        :param v: Initial velocity.
        :param dt: Time step.
    :return: 
        Simulated position and velocity over time.
    """
    # Initialize arrays to store position and velocity
    num_steps = 1000
    positions = np.zeros(num_steps)
    velocities = np.zeros(num_steps)
    positions[0] = x
    velocities [0] = v
    for n in range(1, num_steps):
        # Calculate acceleration from force terms
        a = -(k / m) * positions[n-1] - (gamma / m) * velocities [n-1]
        # Update position
        positions[n] = positions[n-1] + velocities[n-1] * dt + 0.5 * a * dt**2
        # Calculate new acceleration
        a_new = -(k / m) * positions[n] - (gamma / m) * velocities[n-1]
        # Update velocity
        velocities[n] = velocities[n-1] + 0.5 * (a + a_new) * dt
    return positions, velocities
    
# Simulation parameters
mass = 1.0
stiffness = 10.0
damping = 0.5
initial_position = 1.0
initial_velocity = 0.0
time_step = 0.1 # Run simulation

pos, vel = harmonic_oscillator(mass, stiffness, damping, initial_position, initial_velocity, time_step) 
print("Final position:", pos[-1])
print("Final velocity:", vel[-1])

Final position: 1.4012065282144866e-11
Final velocity: 7.488868522582272e-12
In [ ]: