\documentclass[12pt,a4paper]{report}

\usepackage[italian]{babel}
\usepackage{helvet}
\renewcommand{\familydefault}{\sfdefault}
\usepackage{graphics}
 \usepackage{amsfonts}

 \begin{document}
 \parindent=15pt

Consideriamo un sistema termodinamico ed etichettiamo i suoi possibili microstati con l'indice $i$ (\`e facile tradurre i ragionamenti che seguono nel caso di indici continui). Il punto etichettato con $i$ \`e un punto nello spazio delle fasi del sistema se lo descriviamo con la Meccanica Classica e un vettore (a meno di normalizzazioni) in un opportuno spazio di Hilbert se lo descriviamo con la Meccanica Quantistica.

 Consideriamo un numero $n\to \infty$ di copie del sistema e sia $n_i$ il numero di copie che si trovano nello stato $i$.  L'evoluzione temporale del numero di occupazione $n_i$ la possiamo scrivere come una derivata $dn_i/dt$ che sar\`a pari al numero di copie che arrivano nello stato $i$ nell'unit\`a di tempo meno il numero di quelle che lo lasciano. A loro volta, gli arrivi possiamo scriverli come il prodotto della una probabilit\`a nell'unit\`a di tempo $W_{ij}$ che una singola copia arrivi allo stato $i$ proveniente  dallo stato $j$ moltiplicata per il numero di copie presenti in $j$; le partenze saranno analogamente scrivibili come  $W_{ji} n_i$. La matrice ${\bf W}$ non \`e in generale simmetrica.

Si ottiene da questo la {\sl Master Equation}

\begin{equation}
\label{master}
dn_i/dt =\sum_j W_{ij} n_j - \sum_j W_{ji} n_i
\end{equation}

\noindent che descrive l'evoluzione temporale delle occupazioni in funzione del tempo. Per quanto assai plausibile, questa equazione non ha ancora una deduzione rigorosa a partire dalle equazioni della meccanica. Nondimeno, ha avuto numerosissime conferme empiriche e la considereremo valida.

La Master Equation nel caso in cui le probabilit\`a $p_i=n_i/n$ non cambiano, quindi di un sistema in equilibrio termodinamico, diventa

\begin{displaymath}
\sum_j W_{ij} n_j = \sum_j W_{ji} n_i.
\end{displaymath}

Facendo adesso l'ipotesi, detta del {\sl bilancio dettagliato}, che in moltissimi casi pu\`o essere dimostrata, secondo cui non solo le partenze e gli arrivi si bilanciano globalmente, ma che si bilanciano separatamente gli scambi tra ogni coppia di stati $i$ e $j$, l'equazione si semplifica in 

\begin{displaymath}
W_{ij} n_j = W_{ji} n_i\quad \forall i,j .
\end{displaymath}

Il rapporto delle probabilit\`a $W$ \`e quindi pari al rapporto dei numeri di occupazione:

\begin{equation}
\label{rapporto}
W_{ij}/W_{ji}= n_i/n_j.
\end{equation}

Il rapporto a secondo membro dipende dal tipo di ensemble considerato. Per esempio, in un microcanonico il postulato di equiprobabilit\`a a priori d\`a $n_i/n_j=1$ e pertanto la matrice ${\bf W}$ risulta simmetrica: la probabilit\`a di andare da uno stato all'altro \`e identica a quella di tornare indietro. In un canonico, invece, $n_i=n p_i = n e^{-E_i/k_BT}/Z$, dove $Z$ \`e la funzione di partizione, per cui

\begin{equation}
\label{canonico}
W_{ij}/W_{ji}= e^{\left( E_j-E_i\right)/k_BT} = e^{-\delta_{ij}/k_BT}.
\end{equation}

Questa espressione suggerisce un algoritmo per il calcolo della distribuzione $n_i$, detto {\sl algoritmo Montecarlo di Metropolis}:

\begin{enumerate}
\item si prende uno stato $j$;
\item si modifica lo stato in modo casuale ottenendo uno stato tentativo $i$;
\item si calcola la differenza di energie tra gli stati, $\delta_{ij} = E_i-E_j$;
\item se $ \delta_{ij}<0$, ovvero se l'energia cala,  si sostituisce sempre $j$ con $i$ (ovvero "si accetta la mossa");
\item se $ \delta_{ij}>0$ si sostituisce $j$ con $i$ con probabilit\`a $e^{-\delta_{ij}/k_BT}$, altrimenti si lascia lo stato in $j$;
\item si aggiornano i dati di output e si torna al passo 2.
\end{enumerate}

In questa maniera si generano delle sequenze di stati in cui ciascuno dipende solo dal quello immediatamente precedente e non dalla storia passata (senza cio\`e avere memoria), sequenze che quindi costituiscono {\sl catene markoviane}. Queste sequenze non hanno connessione con la reale evoluzione temporale del sistema, se non nel fatto che generano, nel limite di catene infinite e se vale l'ipotesi ergodica, la stessa distribuzione dell'ensemble statistico voluto.

I passi che distinguono il Metropolis dagli altri algoritmi Montecarlo sono i 4 e 5, e possono essere sostituiti con qualunque algoritmo che rispetti la condizione \ref{canonico} (il Metropolis pare essere per\`o in generale il pi\`u efficiente).  Come si vede immediatamente, Metropolis corrisponde a scegliere $W_{ji}=1,\ W_{ij}= e^{-\delta_{ij}/k_BT}$ se  $ \delta_{ij}>0$ e viceversa se $ \delta_{ij}<0$ e la \ref{canonico} \`e rispettata.

Il passo 5 si implementa facilmente tramite il {\sl rejection method}. Si genera un numero casuale $r$ tra 0 e 1 e si confronta con $e^{-\delta_{ij}/k_BT}$; se $r<e^{-\delta_{ij}/k_BT}$ si accetta la mossa.

Il passo 2 pu\`o variare molto a seconda del tipo di sistema: nel caso di sistemi atomici classici per esempio di solito consiste nello spostare uno o pi\`u atomi di una distanza prefissata, o casuale ma di media prefissata, in direzioni prese a caso. Maggiore \`e questa distanza minore sar\`a la probabilit\`a di accettazione, ma meno correlati saranno i dati ottenuti e quindi la statistica sar\`a migliore. Tenendo conto di queste richieste contrastanti, l'algoritmo risulta in generale avere la massima efficienza se le mosse effettuate portano a un tasso di accettazione vicino al 50\%.

{\tt \begin{verbatim}
import numpy as np
import matplotlib.pyplot as plt

def metropolis_harmonic_oscillator(n_steps=100000, delta=0.5, beta=1.0, k=1.0):
    """
    Implementazione dell'algoritmo di Metropolis per un oscillatore armonico.
    
    Parametri:
    n_steps : Numero di iterazioni
    delta   : Ampiezza massima del passo (step size)
    beta    : Inversa della temperatura termodinamica (1/kT)
    k       : Costante elastica della molla
    """
    
    # Inizializzazione
    x_current = 0.0  # Posizione iniziale
    samples = []
    accepted = 0

    for _ in range(n_steps):
        # 1. Proposta di un nuovo stato (distribuzione uniforme attorno a x_current)
        x_proposal = x_current + np.random.uniform(-delta, delta)
        
        # 2. Calcolo della variazione di energia dE = V(x_new) - V(x_old)
        energy_current = 0.5 * k * x_current**2
        energy_proposal = 0.5 * k * x_proposal**2
        diff_energy = energy_proposal - energy_current
        
        # 3. Criterio di accettazione di Metropolis
        if diff_energy <= 0:
            # Se l'energia diminuisce, accetta sempre
            x_current = x_proposal
            accepted += 1
        else:
            # Se l'energia aumenta, accetta con probabilità e^(-beta * dE)
            probability = np.exp(-beta * diff_energy)
            if np.random.rand() < probability:
                x_current = x_proposal
                accepted += 1
        
        samples.append(x_current)

    acceptance_rate = accepted / n_steps
    return np.array(samples), acceptance_rate

# --- Esecuzione e Visualizzazione ---

# Parametri simulazione
beta_val = 2.0
k_val = 1.0
n_samples = 100000

samples, rate = metropolis_harmonic_oscillator(n_steps=n_samples, beta=beta_val, k=k_val)

print(f"Tasso di accettazione: {rate:.2%}")

# Grafico dei risultati
plt.figure(figsize=(10, 6))

# Istogramma dei campioni
count, bins, _ = plt.hist(samples, bins=50, density=True, alpha=0.6, color='skyblue', label='MC Samples')

# Distribuzione teorica di Boltzmann: P(x) = sqrt(beta*k / 2pi) * exp(-0.5 * beta * k * x^2)
x_range = np.linspace(bins[0], bins[-1], 200)
theoretical = np.sqrt(beta_val * k_val / (2 * np.pi)) * np.exp(-0.5 * beta_val * k_val * x_range**2)

plt.plot(x_range, theoretical, 'r-', lw=2, label='Teoria (Boltzmann)')
plt.title(f"Metropolis Monte Carlo: Oscillatore Armonico (beta={beta_val})")
plt.xlabel("Posizione x")
plt.ylabel("Densità di probabilità")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
\end{verbatim}}

\end{document}