relazione_sistemi_complessi.tex

\documentclass[10pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[italian]{babel}

%\usepackage[
%  margin=3cm,
%  includefoot,
%  footskip=30pt,
%]{geometry}

\usepackage{listings}
\lstset{basicstyle=\small\ttfamily}

\usepackage{amsmath}
\usepackage{mathtools}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{siunitx}
\usepackage{physics}
\usepackage{dsfont}

\mathtoolsset{showonlyrefs,showmanualtags}

\usepackage{graphicx}
\usepackage{subfig}
\usepackage{wrapfig}
\usepackage{sidecap}
\usepackage{booktabs}
\usepackage{hyperref}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[theorem]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{defn}[theorem]{Definizione}

\newtheorem{prop}{Proposizione}[section]
\newtheorem{dimo}{Dimostrazione}[section]

%%% BackEnd Bibliografico
%\usepackage{textcomp}
%\usepackage[autostyle]{csquotes}
%\usepackage[
%        backend=biber,
%        %bibstyle=numeric,
%        %sorting=ynt
%    ]{biblatex}
%\addbibresource{bibliografia.bib}
%\nocite{*}
%%%

%%%%%% TESTO EFFETTIVO

\title{Relazione di Fisica dei Sistemi Complessi}
\author{Carlo Emilio Montanari}

\begin{document}

\maketitle

\tableofcontents

\section{Concetti preliminari}
In questa sezione si trattano brevemente esempi e proprietà matematiche degli oggetti in uso, necessari alla trattazione del problema in questione di un oscillatore armonico con forzante stocastica.

\subsection{Derivazione dell'equazione di Fokker-Planck}
Consideriamo l'equazione di Langevin nella sua forma più generale:
\begin{equation}
	\dot{x} = a(x,t) + b(x,t)\xi(t)
	\label{eq:stoc_base}
\end{equation}
abbiamo che per una assegnata realizzazione del rumore bianco \(\xi(t)\) (a media nulla e varianza unitaria) abbiamo la corrispondente equazione di Liouville stocastica per \(\rho(x,t;\xi)\) nella forma
\begin{equation}
	\pdv{\rho}{t} + \pdv{x}(a\rho) + \xi(t)\pdv{x}(b\rho)=0
\end{equation}
Riscrivendola nella forma compatta
\begin{equation}
	\pdv{\rho}{t} = \vb{A}\rho +\xi(t)\vb{B}\rho
\end{equation}
avendo definito gli operatori
\begin{equation}
	\vb{A} = -\pdv{}{x} a \quad\quad \vb{B} = - \pdv{}{x} b 
\end{equation}
abbiamo la seguente proposizione:
\begin{prop}
	Il valore medio \(\langle\rho\rangle\) rispetto a tutte le possibili realizzazioni \(\xi\) di \(\rho(x,t;\xi)\) è una densità di probabilità che soddisfa l'equazione di Fokker-Planck:
	\begin{equation}
		\pdv{}{t}\langle\rho\rangle = \vb{A}\langle\rho\rangle + \frac{1}{2}\vb{B}^2\langle\rho\rangle
	\end{equation}
\end{prop}

\paragraph{Dimostrazione}
Consideriamo lo sviluppo di Dyson per una equazione operatoriale nella forma
\begin{equation}
	\pdv{\vb{U}}{t} = \vb{H}(t)\vb{U} \quad\quad\quad \vb{U}(0) = I
\end{equation}
tale sviluppo si ottiene a partire dalla forma integrale
\begin{equation}
	\vb{U}(t) = I + \int_0^t \vb{H}(s)\vb{U}(s)\,ds
\end{equation}
Questa forma integrale può essere iterata nella seguente forma
\begin{align}
	\vb{U}(t) &= I + \sum_{n=1}^{\infty}\int_0^t dt_1 \int_0^{t_1} dt_2 \cdots \int_0^{t_{n-1}}dt_n \vb{H}(t_1) \vb{H}(t_2)\cdots \vb{H}(t_n) = \\
	&= I + \sum_{n=1}^{\infty} \frac{1}{n!} \int_0^t dt_1 \int_0^{t} dt_2 \cdots \int_0^{t}dt_n \mathds{T}(\vb{H}(t_1) \vb{H}(t_2)\cdots \vb{H}(t_n))\\
	&= \exp(\mathds{T} \int_0^y \vb{H}(s)\,ds)
\end{align}
Dove \(\mathds{T}\) è l'operatore di ordinamento temporale. Per giustificare l'ultima uguaglianza basta considerare che (caso per \(n=2\))
\begin{multline}
	\frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \mathds{T}(\vb{H}(t_1) \vb{H}(t_2)) \\
	=\frac{1}{2} \int_0^t dt_1 \int_0^{t_1} dt_2 \vb{H}(t_1) \vb{H}(t_2) + \frac{1}{2} \int_0^t dt_1 \int_{t_1}^{t} dt_2 \vb{H}(t_2) \vb{H}(t_1)\\
	=\int_0^t dt_1 \int_0^{t_1} dt_2 \vb{H}(t_1) \vb{H}(t_2)
\end{multline}
dove l'ultima uguaglianza è stata ottenuta cambiando l'ordine di integrazione del secondo termine nella seconda riga e ridefinendo le variabili mute. Questo procedimento è applicabile a tutti i valori superiori di \(n\).

Consideriamo ora l'operatore di evoluzione \(\vb{U}_{A}\), generato da \(\vb{A}\)
\begin{equation}
	\pdv{\vb{U}_A}{t} = \vb{A}\vb{U}_A \quad\quad\quad \vb{U}_A(0) = \vb{I}
\end{equation}
Ponendo \(\vb{U} = \vb{U}_A\) si ha come conseguente risultato che
\begin{equation}
	\vb{U}_{A} \pdv{\vb{U}_I}{t} = \xi \vb{B}\vb{U}_A \vb{U}_I
\end{equation}
che, ponendo la notazione \(\vb{B}_I = \vb{U}_A^{-1} \vb{B}\vb{U}_A\) si scrive
\begin{equation}
	\pdv{\vb{U}_I}{t} = \xi \vb{B}_I \vb{U}_I
\end{equation}
Da qui, vale la soluzione per \(\vb{U}_I\)
\begin{equation}
	\vb{U}_I(t) = \exp\left(\mathds{T}\int_0^t \vb{B}_I(s)\xi(s)\right)
\end{equation}
Si può quindi calcolare la media della soluzione di \(\vb{U}_I\) rispetto alle realizzazioni di \(\xi \) usando lo sviluppo in serie della soluzione e considerando la forma del valor medio dei prodotti del rumore bianco, che è nullo quando il numero di fattori è dispari. Si ha quindi
\begin{align}
	\langle\vb{U}_I\rangle &= \sum_{n=0}^\infty \frac{(2n-1)!!}{(2n)!}\int_0^t dt_1 \cdots \int_0^t dt_{2n} \mathds{T}(\vb{B}_I(t_1)\cdots\vb{B}_I(t_{2n})) \delta(t_1 + t_{n+1})\delta(t_2 + t_{n+2})\cdots\delta(t_n + t_{2n}) \\
	&= \sum_{n=0}^\infty \frac{1}{2^n n!} \mathds{T}\left(\int_0^t dt_1 \cdots \int_0^t dt_n \vb{B}_I^2(t_1)\cdots\vb{B}_I^2(t_n)\right) \\
	&= \exp(\frac{1}{2}\mathds{T}\int_0^t ds\, \vb{B}_I^2(s))
\end{align}
Se ne ricava quindi che
\begin{equation}
	\pdv{\langle\vb{U}_I\rangle}{t} = \frac{1}{2} \vb{B}_I^2\langle\vb{U}_I\rangle = \frac{1}{2} \vb{U}_A^{-1}\vb{B}^2\vb{U}_A\langle\vb{U}_I\rangle
\end{equation}
Avendo infine che \(\langle\vb{U}\rangle = \vb{U}_A\langle\vb{U}_I\rangle \) si ha che
\begin{equation}
	\pdv{\langle \vb{U}\rangle}{t} = \pdv{\vb{U}}{t}\langle\vb{U}_I\rangle + \vb{U}_A \pdv{\langle\vb{U}_I\rangle}{t} = A\vb{U}_A\langle\vb{U}_I\rangle + \frac{1}{2} \vb{B}^2_I\vb{U}_A\langle\vb{U}_I\rangle = A\langle\vb{U}_A\vb{U}_I\rangle + \frac{1}{2}\vb{B}^2\langle\vb{U}_A\vb{U}_I\rangle
\end{equation}
che dimostra il risultato.



\subsection{Diffusione su superficie toroidale}
\label{sec:diffusione_toro}
Definiamo per prima cosa la seguente notazione per il rumore di Wiener ed il suo integrale (ove non indicato diversamente, si fa uso sempre della definizione di integrale secondo Ito):
\begin{equation}
	w(t) = \int_0^t \xi(s)\,ds \quad\quad w_1(t)=\int_0^t w(t)\,ds
\end{equation}
Abbiamo che possiamo riscrivere \(w_1(t)\) nella forma
\begin{equation}
	w_1(t) = \int_0^t \int_0^s \xi(u)\,du\,ds = \int_0^t \xi(u)\int_u^s\,ds\,du = \int_0^t (t-u)\xi(u)\,du
\end{equation}
Rappresentando quindi i processi \(w(t)\) e \(w_1(t)\) nella forma standard \(\int_0^t K(t,s)\xi(s)\,ds\), si ha che \(K=1\) per \(w(t)\) e che \(K=t-s\) per \(w_1(t)\). Dalla teoria abbiamo che \(\sigma^2 = t\) per il rumore di Wiener \(w(t)\) mentre, per il suo integrale \(w_1(t)\), abbiamo che:
\begin{equation}
	\sigma^2(t) = \int_0^t (t-u)^2\,du = \frac{t^3}{3}.
\end{equation}

Consideriamo a questo punto l'esempio più semplice di diffusione angolare per un sistema anisocrono. Sia \(x\) una variabile angolare divisa per \(2\pi\) definita sul toro \([0, 1]\). Sia poi
\begin{equation}
	\dot{x} = \omega + y \quad\quad \dot{y} = \epsilon\xi(t)
	\label{eq:langevin_toro}
\end{equation}
Se si ha una soluzione di diffusione per \(x\) definita su \(\mathbf{R}\) nella forma \(\rho(x,t)\) è possibile ottenere la soluzione sul toro \(\rho_T(x,t)\) tramite la periodicizzazione
\begin{equation}
	\rho_T(x,t) = \sum_{k=-\infty}^{+\infty} \rho(x+k,t).
\end{equation}
La soluzione di \eqref{eq:langevin_toro} su \(\mathbf{R}\) è data da
\begin{equation}
	y = y_0 + \epsilon w(t) \quad\quad x = x_0 + (\omega + \epsilon y_0) t + \epsilon w_1(t)
\end{equation}
Indicando quindi la parte deterministica di \(x\) con \(\langle x(t) \rangle = x_0 + (\omega + \epsilon y_0) t\), si può scrivere la densità di probabilità su \(\mathbf{R}\) nella forma
\begin{equation}
	\rho(x,t) = \frac{\exp{\frac{-(x-\langle x(t)\rangle)^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \quad\quad \sigma^2 = \epsilon^2 \frac{t^3}{3}
\end{equation}
da cui consegue che
\begin{equation}
	\rho_T(x,t) = \sum_{n=-\infty}^{+\infty} \frac{\exp{\frac{-(x-\langle x(t)\rangle + n)^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}}.
\end{equation}
Rappresentando quindi la soluzione in forma di serie di Fourier
\begin{equation}
	\rho_T(x,t) = \sum_{k=-\infty}^{+\infty} f_k(t)e^{2\pi i k x}
\end{equation}
con i coefficienti \(f_k\) dati da
\begin{align}
	f_k &= \int_0^1 e^{-2\pi i k x} \rho_T(x,t)\,dx \\
	&= \sum_{n=-\infty}^{+\infty} \frac{1}{\sqrt{2\pi\sigma^2}} \int_0^1 \exp{\frac{-(x-\langle x(t) \rangle + n)^2}{2\sigma^2} -2\pi i k x}\,dx
\end{align}
portando avanti i conti e supponendo i coefficienti \(f_k\) reali, si ottiene
\begin{equation}
	\rho_T(x,t) = 1+2\sum_{k=1}^\infty e^{-2\pi^2\sigma^2k^2} f_k \cos(2\pi k (x-\langle x(t) \rangle)).
\end{equation}
Il punto più importante di questo esempio è notare come \(\rho_T\) rilassa alla distribuzione uniforme come \(e^{-\sigma^2}\), tenendo sempre a mente che, in questo caso,
\begin{equation}
	\sigma^2 \propto t^{3}.
\end{equation}

\subsubsection{Considerazioni per sistemi hamiltoniani}
Consideriamo due casi semplici ma significativi di sistemi hamiltoniani in coordinate azione-angolo perturbati stocasticamente nella forma
\begin{equation}
	H = H_0(J) + \epsilon H_1(\theta, J)\xi(t).
\end{equation}
\paragraph{Esempio 1.}
Consideriamo il caso in cui
\begin{equation}
	H = H_0(J) + \epsilon J \xi(t) 
\end{equation}
abbiamo che
\begin{align}
	\dot{J} &= 0 \\
	\dot{\theta} &= \omega(J) + \epsilon\xi
\end{align}
dove \(\omega(J) = \partial H_0/\partial J\). Abbiamo quindi che \(\theta\) descrive un semplice processo di Wiener la cui soluzione può essere scritta nella forma
\begin{equation}
	\theta = \theta_0 + \omega t + \epsilon w(t)
\end{equation}
da cui consegue una varianza tale che
\begin{equation}
	\sigma_{\theta}^2 \propto t
\end{equation}

\paragraph{Esempio 2.}
Consideriamo il caso in cui
\begin{equation}
	H = H_0(J) - \epsilon \theta \xi(t)
\end{equation}
abbiamo che
\begin{align}
	\dot{J} &= \epsilon\xi \\
	\dot{\theta} &= \omega(J)
\end{align}
da cui consegue immediatamente che
\begin{equation}
	J = J_0 + \epsilon w(t)
\end{equation}
se poi andiamo a fare uno sviluppo di Taylor per \(\dot{\theta}\) abbiamo che
\begin{equation}
	\dot{\theta} = \omega(J_0) + \omega'(J_0)\epsilon w(t) + o(\epsilon^2)
\end{equation}
che, integrato, porta a
\begin{equation}
	\theta = \theta_0 + \omega(J_0)t + \epsilon \omega'(J_0) w_1(t)
\end{equation}
da cui, sulla base delle osservazioni fatte a inizio sezione, consegue una varianza di \(\theta\) tale che
\begin{equation}
	\sigma_{\theta}^2 \propto t^3
\end{equation}
mentre per l'azione \(J\) si ha una varianza
\begin{equation}
	\sigma_{J}^2 \propto t
\end{equation}
Considerato quindi che il tempo di rilassamento della variabile angolare è sempre molto minore del tempo di rilassamento della variabile d'azione, si giustifica quindi per tempi sufficientemente lunghi l'approssimazione di media sulla variabile angolare per le distribuzioni di probabilità.


\subsection{Sistemi hamiltoniani con perturbazione stocastica}
\label{sec:hamiltoniana_stocastica}
Consideriamo un sistema hamiltoniano perturbato, descritto in coordinate azione-angolo con origine in un punto critico stabile di \(H_0\) e le coordinate \((\theta, J)\) definite nella regione delimitata dalla separatrice. Poniamolo nella forma
\begin{equation}
	H = H_0(J) + \epsilon H_1(\theta, J)\xi(t).
\end{equation}
Si ha che la densità di probabilità \(\rho(\theta, J, t)\) soddisfa l'equazione
\begin{equation}
	\pdv{\rho}{t} + \Omega \pdv{\rho}{\theta} = \frac{\epsilon^2 \sigma^2}{2}[H_1, [H_1,\rho]] 
\end{equation}
dove \([,]\) è la parentesi di Poisson e \(\Omega(J) = \partial H_0/\partial J\).

Se si tiene poi conto che, come mostrato nelle sottosezioni precedenti, l'angolo in quanto variabile toroidale rilassa ad una distribuzione uniforme di equilibrio con un tempo \(\tau_{\theta} \propto \epsilon^{-2/3}\) molto piccolo rispetto alla scala dei tempi di rilassamento \(\tau_{J} \propto \epsilon^{-2}\) dell'azione, è possibile semplificare di molto quest'ultima equazione. Si ha che per \(t\gg\tau_{\theta}\) l'equazione per la distribuzione mediata sull'angolo \(\rho(J,t)\) soddisfa l'equazione
\begin{equation}
	\pdv{}{t}\rho(J,t) = \pdv{}{J}D(J)\pdv{}{J}\rho(J,t)
\end{equation}
dove
\begin{equation}
	D(J) = \frac{\epsilon^2 \sigma^2}{2}\frac{1}{2\pi}\int_0^{2\pi}\left(\pdv{H_1}{\theta}\right)^2 \,d\theta.
\end{equation}

\paragraph{Dimostrazione:}
Partiamo dall'equazione di Liouville stocastica:
\begin{equation}
	\pdv{\rho}{t} + \Omega(J) \pdv{\rho}{\theta} + \epsilon\xi(t)[\rho,H_1] = 0
	\label{eq:dimo_liouville}
\end{equation}
Scriviamo ora \(\rho\) nella forma \(\rho = \rho_0 + \epsilon\rho_1\) dove \(\rho_0\) è la componente media e \(\rho_1\) la parte fluttuante a media nulla, considerando anche che \(\langle\xi\rangle = 0\) abbiamo che il valor medio di \eqref{eq:dimo_liouville} è
\begin{equation}
	\pdv{\rho_0}{t} + \Omega(J)\pdv{\rho_0}{\theta} + \epsilon^2[\langle\xi(t)\rho_1\rangle, H_1] = 0
	\label{eq:dimo_media}
\end{equation}
Sottraendo ora \eqref{eq:dimo_media} da \eqref{eq:dimo_liouville}, si ottiene l'equazione
\begin{equation}
	\pdv{\rho_1}{t} + \Omega(J) \pdv{\rho_1}{\theta} = -\xi(t)[\rho_0, H_1] + O(\epsilon)
\end{equation}
che desideriamo risolvere per \(\rho_1\) in modo da poter sostituire in \eqref{eq:dimo_media}. Per riuscirci eseguiamo il cambio di variabili
\begin{align}
	\theta &\to \theta - \Omega \tau \\
	t &\to t - \tau
\end{align}
che ci permette di scrivere
\begin{equation}
	\dv{}{\tau}\rho_1(\theta - \Omega \tau, J, t - \tau) = \xi(t- \tau)[\rho_0, H_1](\theta - \Omega \tau, J, t - \tau)
\end{equation}
che possiamo integrare da \(\tau = 0\) a \(\tau = t\), ottenendo
\begin{equation}
	\rho_1(\theta, J, t) = -\int_0^t[\rho_0, H_1](\theta - \Omega \tau, J, t - \tau)\xi(t-\tau)\,d\tau
\end{equation}
dove abbiamo sfruttato il fatto che \(\rho_1(\theta, J, 0)=0\). Moltiplicando poi entrambi i membri per \(\xi(t)\) e calcolando il valor medio rispetto a tutte le realizzazioni del rumore, abbiamo nel caso di un rumore di Wiener che
\begin{equation}
	\langle \xi(t)\rho_1(\theta, J, t) \rangle = -\frac{1}{2}\sigma^2 [\rho_0, H_1](\theta, J, t)
\end{equation}
Andando quindi a sostituire questo risultato in \eqref{eq:dimo_media}, andiamo ad ottenere, a meno di termini di ordine \(\epsilon^3\), l'equazione per la densità media
\begin{equation}
 	\pdv{\rho_0}{t} + \Omega (J) \pdv{\rho_0}{\theta} = \frac{\epsilon^2 \sigma^2}{2}[[\rho_0, H_0], H_0].
\end{equation} 
Come abbiamo visto nella sezione \ref{sec:diffusione_toro}, per tempi sufficientemente lunghi abbiamo che in buona approssimazione \(\rho = \rho(J,t)\) non dipende dall'angolo e fa avere come doppia parentesi di Poisson
\begin{equation}
	[[\rho_0, H_0], H_0] = \pdv{}{J}\left[\left(\pdv{H_1}{\theta}\right)^2 \pdv{\rho_0}{J}\right] - \pdv{}{\theta}\left[\pdv{\rho_0}{J}\pdv{H_1}{\theta}\pdv{H_1}{J}\right]
\end{equation}
Prendendo quindi la media angolare di quest'ultima equazione, il secondo termine del lato destro si annulla e si trova l'espressione del coefficiente di diffusione \(D(J)\) da dover integrare lungo l'intero angolo.

\subsection{Fluttuazione e Dissipazione}
\label{sec:fluttuazione_dissipazione}
Consideriamo ora la seguente equazione di Langevin
\begin{align}
	\dv{x}{t} &= \frac{p}{m}\\
	\dv{p}{t} &= F_{ex} - \gamma p + \epsilon\xi(t)
\end{align}
Abbiamo che la funzione di densità nello spazio delle fasi \(\rho(x,p,t)\) soddisfa l'equazione di Fokker-Planck
\begin{equation}
	\pdv{\rho}{t} + \pdv{x} (\frac{p}{m}\rho) + \pdv{p} (F_{ex}\rho - \gamma p \rho) = \epsilon^2\pdv[2]{\rho}{p}
\end{equation}
Se la forza esterna \(F_{ex}\) è conservativa ed ammette un potenziale \(V\), è possibile trovare una soluzione analitica di equilibrio per questa equazione di Fokker-Planck. La soluzione di equilibrio si scrive nella forma:
\begin{equation}
	\rho(x,p) = Z^{-1} \exp{-\frac{H}{k_B T}}
\end{equation}
dove si è posto
\begin{equation}
	H = \frac{p^2}{2m} + V(x) \quad\quad k_b T = \frac{\epsilon^2}{2 m \gamma}
\end{equation}
e dove con \(Z\) si indica la costante di normalizzazione tale che \(\int f\,dx\,dp = 1\). È evidente che la soluzione di equilibrio corrisponde alla distribuzione di Maxwell Boltzmann.

\subsection{Averaging Principle}
\label{subsec:averaging_principle}
Se una piccola perturbazione (e.g.\ un processo stocastico con \(\epsilon \) piccolo) viene applicata ad un sistema conservativo, si ha che le quantità che prima erano integrali primi del moto cominciano ad evolvere lentamente. Su scale di tempo di ordine 1 tali evoluzioni sono piccole, su tempi di ordine \(1/\epsilon \) l'evoluzione può essere non trascurabile.

Tramite l'\textit{averaging principle}, è possibile scrivere una equazione di evoluzione contenente solamente tali variabili lente, trascurando le variabili veloci. Questo, va ribadito, è un principio fisico e non un teorema e, in quanto tale, non ha una formulazione matematica rigorosa.

Assumiamo un sistema in coordinate azione-angolo in cui le coordinate angolari appartengono ad una superficie toroidale. Assumiamo inoltre che le equazioni del moto non perturbato possano essere scritte nella forma:
\begin{align}
	\dot{I} &= 0,\\
	\dot{\theta} &= \omega(I)
\end{align}
Una piccola perturbazione al sistema fa sorgere nuovi termini nelle equazioni del moto, questi termini assumono la forma:
\begin{align}
	\dot{I} &= \epsilon f(I, \theta, \epsilon),\\
	\dot{\theta} &= \omega(I) + \epsilon g(I, \theta, \epsilon)
\end{align} 
dove $f$ e $g$ hanno periodo $2\pi$ in $\theta$ ed $\epsilon$ è piccolo. La variabile azione $I$ viene detta \textit{variabile lenta} mentre la variabile angolo $\theta$ viene detta \textit{variabile veloce}.

Essendo di norma interessante conoscere il comportamento nel tempo della variabile lenta, l'averaging principle consiste nel rimpiazzare il sistema di equazioni perturbate con il sistema mediato
\begin{equation}
	\dot{J} = \epsilon F(J), \qquad F(J) = \frac{1}{2\pi} \oint_{\mathds{T}} f(J,\theta,0)\,d\theta
\end{equation}
per descrivere l'evoluzione approssimata delle variabili lente del sistema lungo tempi di ordine $1/\epsilon$. 

\subsection{Invarianza adiabatica nell'oscillatore armonico}
Sia \(H(q,p;\lambda)\) una hamiltoniana con un singolo grado di libertà. Supponiamo che \(\lambda(t)\) vari lentamente nel tempo, dove per lentamente si intende matematicamente che esiste un \(\epsilon\) tale che
\begin{equation}
	\frac{1}{n!}\abs{\dv[n]{\lambda}{t}} \leq \epsilon^n
\end{equation}
Possiamo per esempio lavorare semplicemente con \(\lambda = \epsilon t\) con \(\epsilon\) molto piccolo.

Si definisce un invariante adiabatico una funzione \(J(p,q;\lambda)\) del sistema \(H\) se questa è tale che, per ogni \(\epsilon > 0\), esiste un \(\epsilon_0>0\) tale che, se \(\epsilon < \epsilon_0\), \(0<t<(1/\epsilon)\) si ha che
\begin{equation}
	\abs{J(q(t),p(t);\epsilon t) - J(q(0), p(0);0)}<c\epsilon^\alpha \quad\quad (\alpha > 0)
\end{equation}
Per esempio, se abbiamo un invariante adiabatico \(J\), abbiamo, per un tempo \(t<(1/\epsilon)\), \(\alpha=1\), che \(J(t)-J(0) = \mathcal{O}(\epsilon)\).

È possibile dimostrare che ogni sistema unidimensionale ammette un invariante adiabatico e che la variabile di Azione risulta essere un invariante adiabatico. Ricordiamo che quando ci si riferisce a variabili Azione-Angolo in un sistema hamiltoniano \(H(q,p)\) unidimensionale, ci si riferisce ai seguenti integrali
\begin{align}
	J &= \frac{1}{2\pi} \oint_{H=E} p\,dq \\
	\theta &= \pdv{}{J} \int_{H=E} p\,dq
\end{align}
dove la curva \(\gamma(E)\) è una curva di livello chiusa nello spazio delle fasi all'interno di una regione delimitata da una separatrice. Ricordiamo inoltre che il cambio di coordinate \((p,q)\to(J,\theta)\) è canonico e che possiamo definire la frequenza \(\omega\) nella forma \(\omega = \partial H / \partial J\).

Analizziamo quindi il caso di un oscillatore armonico unidimensionale con frequenza \(\omega(\lambda)\) lentamente variante nel tempo. Consideriamo \(\lambda = \epsilon t\) e \(\omega(\lambda) = \lambda(\omega_1 - \omega_0) + \omega_0\).\\
Prima di considerare la perturbazione, ricordiamo che l'hamiltoniana dell'oscillatore armonico è
\begin{equation}
	H = \frac{p^2}{2} + \omega^2\frac{q^2}{2}
\end{equation}
e che, per energia \(E\) fissata, il moto descrive una curva di livello di forma ellittica che rende immediato il calcolo della variabile Azione
\begin{equation}
	J = \frac{E}{\omega} = \frac{1}{2\omega}(p^2 + \omega^2 q^2)
\end{equation}
Introducendo quindi la perturbazione causata da \(\lambda\) (consideriamo \(\omega_0 = 1\) e \(\omega_1 = 2\), avendo quindi \(\omega = 1+\lambda\)), abbiamo
\begin{equation}
	J = \frac{q^2(1+\epsilon t)}{2} + \frac{p^2}{2(1+\epsilon t)}
\end{equation}
che sviluppando con Taylor diventa
\begin{equation}
	J = \frac{q^2}{2} + \frac{p^2}{2} + \epsilon\left(\frac{q^2 t}{2}-\frac{p^2 t}{2}\right) + \mathcal{O}(\epsilon^2)
\end{equation}
Per poter dire che \(J\) è un invariante adiabatico, bisogna dimostrare che, lungo un periodo,
\begin{equation}
	\Delta J = \mathcal{O}(\epsilon^2)
\end{equation}
Considerando il periodo di una oscillazione completa \(T\),
\begin{equation}
	\Delta J = J(T) - J(0) = \epsilon\left(\frac{q^2}{2} - \frac{p^2}{2}\right)T + \mathcal{O}(\epsilon^2)
\end{equation}
che possiamo riesprimere nella forma
\begin{equation}
	\Delta J = \epsilon(H(T) - H(0))T - \epsilon^2 T^2 \frac{q^2}{2} + \mathcal{O}(\epsilon^2)
\end{equation}
e, avendo che in un periodo \(H(T) - H(0) = \epsilon T \frac{q^2}{2}\), si dimostra che
\begin{equation}
	\Delta J = \mathcal{O}(\epsilon^2)
\end{equation}
e si ha quindi che \(J\) è un invariante adiabatico.


\subsection{Integrazione numerica di una equazione di Langevin}
\label{sec:integratore}
Consideriamo sempre una equazione di Langevin nella seguente forma:
\begin{align}
	dq(t) &= p(t)dt \\
	dp(t) &= f(q(t))dt - \gamma p(t)dt + \sigma dw(t)
\end{align}
dove \(w(t)\) è un processo di Weiner. Vogliamo implementare una discretizzazione valida di questa equazione, riuscendo a mantenere le proprietà statistiche principali causate dalla presenza del processo di Weiner e che mantenga una accuratezza di ordine \(\Delta t^2\). Se cominciamo integrando il sistema da \(t\) a \(t + \Delta t\)
\begin{align}
	q(t+\Delta t) &= q(t) + \int_t^{t+\Delta t} p(s)\,ds\\
	p(t+\Delta t) &= p(t) + \int_t^{t+\Delta t} f(q(s))\,ds -\gamma \int_t^{t+\Delta t} p(s)\,ds + \,\, \sigma[w(t+\Delta t)-w(t)]
\end{align}

Entrambe le equazioni contengono integrali in \(p\) e in \(f\) che producono termini di ordine \(\Delta t\) e superiore nella discretizzazione. Cominciamo ad analizzare l'equazione per \(p\), che può essere usata all'ordine più basso:
\begin{equation}
	p(s) \approx p(t) + (s-t)f(q(t)) -(s-t)\gamma p(t) + \sigma[w(s)-w(t)]
\end{equation}
dove \(t\leq s\leq t+\Delta t\). Integrando l'equazione in questa forma da \(t\) a \(t+\Delta t\) otteniamo:
\begin{equation}
	\int_t^{t+\Delta t} p(s)\,ds = \Delta t p(t) + C(t)
\end{equation}
dove si è definito
\begin{equation}
	C(t) \equiv \frac{\Delta t^2}{2}[f(q(t))-\gamma p(t)] + \sigma \int_t^{t+\Delta t}[w(s)-w(t)]\,ds
	\label{eq:c_t}
\end{equation}
questo può essere usato come approssimazione degli integrali su \(p\) che compaiono nelle equazioni precedenti.

Per approssimare l'integrale su \(f\), cominciamo sfruttando la relazione
\begin{equation}
	\dv{f}{t} = \pdv{f}{q} \dot{q} = \pdv{f}{q} p
\end{equation}
ottenendo quindi
\begin{equation}
	f(q(s)) = f(q(t))+\int_t^s p(k)f'(q(k))\,dk \approx f(q(t))+(s-t)p(t)f'(q(t))
\end{equation}
e, come prima, integriamo questa da \(t\) a \(t+\Delta t\), ottenendo
\begin{equation}
	\int_t^{t+\Delta t} f(q(s))\,ds = \Delta t f(q(t)) + \frac{\Delta t^2}{2}p(t)f'(q(t)) = \Delta t \frac{f(q(t+\Delta t))+ f(r(t))}{2}
\end{equation}
dove l'ultima uguaglianza è stata ottenuta tramite la seguente relazione:
\begin{equation}
	f(q(t+\Delta t)) = f(q(t)) + f'(q(t))(q(t+\Delta t)-q(t)) = f(q(t)) + f'(q(t))p(t)\Delta t.
\end{equation}

Per computare correttamente il termine \(w(t+\Delta t) - w(t)\), sfruttiamo la relazione
\begin{equation}
	w(t+\Delta t) - w(t) = \sqrt{\Delta t}\,\xi \label{eq:xi_wiener}
\end{equation}
resta però da calcolare l'integrale della differenza tra due funzioni di Wiener nell'espressione \eqref{eq:c_t} di \(C(t)\). Sfruttando il fatto che
\begin{equation}
	\langle w(s)w(s') \rangle = \int_0^s \int_0^{s'} \langle \eta(t)\eta(t') \rangle\,dt\,dt' = \int_0^s \int_0^{s'} \delta(t'-t) \,dt\,dt' = \min(s,s')
\end{equation}
si può dimostrare che:
\begin{align}
	\left\langle (w(t+\Delta t)-w(t))\int_t^{t+\Delta t} (w(s')-w(t))\,ds' \right\rangle &= \frac{\Delta t^2}{2} \label{eq:cosaaaaaa}\\
	\left\langle \int_t^{t+\Delta t} (w(s)-w(t))\,ds \int_t^{t+\Delta t} (w(s')-w(t))\,ds' \right\rangle &= \frac{\Delta t^3}{3} \label{eq:da_rispettare}
\end{align}
Possiamo quindi fare in modo che l'equazione \eqref{eq:da_rispettare} sia rispettata esprimendo
\begin{equation}
	\int_t^{t+\Delta t} (w(s) - w(t))\,ds = \Delta t^{3/2}\eta
	\label{eq:forma_quasi_finale}
\end{equation}
con \(\eta\) numero random tale che \(\langle\eta\rangle = 0\) e \(\langle\eta^2\rangle = 1/3\). Sfruttando quindi \eqref{eq:xi_wiener} e \eqref{eq:da_rispettare} in \eqref{eq:forma_quasi_finale}, otteniamo \(\langle \xi \eta \rangle = 1/2\), ossia, i due numeri random sono inevitabilmente correlati. Per risolvere questo problema riformuliamo \(\eta\) in due variabili random non correlate:
\begin{equation}
	\eta = a \xi + b \theta
\end{equation}
tali che \(\langle\theta\rangle = 0\), \(\langle\theta^2\rangle = 1\) e \(\langle\xi\theta\rangle = 0\). Determiniamo quindi \(a\) e \(b\) imponendo che \(\langle \xi \eta\rangle = 1/2\) e che \(\langle\eta^2\rangle = 1/3\) per ottenere quindi \(a = 1/2\) e \(b = 1/(2\sqrt{3})\). Otteniamo finalmente:
\begin{equation}
	\int_t^{t+\Delta t} (w(s) - w(t))\,ds = \Delta t^{3/2} \left( \frac{\xi}{2} + \frac{\omega}{2\sqrt{3}} \right)
\end{equation}
Mettendo insieme tutti i risultati raccolti fino ad ora, si ottiene \textbf{L'integratore al secondo ordine dell'equazione di Langevin:}
\begin{align}
 	C(t) &= \frac{\Delta t^2}{2} [f(q(t)) - \gamma p(t)] + \sigma \Delta t^{3/2} \left(\frac{\xi(t)}{2} + \frac{\theta(t)}{2\sqrt{3}} \right) \\
 	q(t+\Delta t) &= q(t) + \Delta t p(t) + C(t) \\
 	p(t+\Delta t) &= p(t) + \frac{\Delta t}{2}[f(q(t+\Delta t)) + f(q(t))] - \Delta t \gamma p(t) + \sigma \sqrt{\Delta t}\,\xi(t) - \gamma C(t)
\end{align}
Dove \(\xi\) e \(\theta\), nel contesto di una simulazione in Python, vengono espressi ad ogni step temporale estraendo un valore casuale da una distribuzione normale di media nulla e varianza unitaria tramite la funzione \lstinline{numpy.random.normal()}.

\subsection{Integrazione di una equazione di Fokker-Planck con metodo di Crank-Nicolson}

In questa sezione si tratterà in breve l'implementazione di un integratore di Crank-Nicolson per una equazione differenziale nella forma:
\begin{equation}
	\pdv{u}{t} = - A \pdv{u}{x} + B \pdv[2]{u}{x}
\end{equation}
dove consideriamo i termini $A$ e $B$ come valori costanti e condizioni al contorno naturali (i.e. barriere assorbenti).

Per convertire ogni termine dell'equazione in forma di differenziale finito scriviamo:
\begin{align}
	\pdv{u}{t} &\Rightarrow u_t \Rightarrow \frac{u_m^{n+1}-u_m^n}{k}\\
	\pdv[2]{u}{x} &\Rightarrow u_{xx} \Rightarrow \frac{(u_{m+1}^{n+1}-2u_m^{n+1}+u_{m-1}^{n+1}) +(u_{m+1}^{n}-2u_m^{n}+u_{m-1}^{n})}{2(h)^2}\\
	\pdv{u}{x} &\Rightarrow u_x \Rightarrow \frac{(u_{m+1}^{n+1}-u_{m-1}^{n+1})+(u_{m+1}^n -u_{m-1}^n)}{4(h)}
\end{align}
definendo poi due costanti per semplificare i conti
\begin{equation}
	\alpha = \frac{Ak}{4(h)}, \quad \beta = \frac{Bk}{2 h^2}
\end{equation}
otteniamo
\begin{multline}
	(\alpha - \beta)(u_{m+1}^{n+1}) + (1 + 2\beta)(u_m^{n+1}) + (-\alpha -\beta)(u_{m-1}^{n+1}) \\= (-\alpha + \beta)(u_{m+1}^n) +(1 -2\beta)(u_m^n) + (\alpha + \beta)(u_{m-1}^n)
\end{multline}
con a sinistra dell'uguale i termini non noti ed a destra dell'uguale i termini noti. È evidente come il problema può venir posto in forma vettoriale
\begin{equation}
	\left(L\right)\vb{u}^{n+1} = \left(R\right)\vb{u}^n + \vb{b},
\end{equation}
dove
\begin{align}
	L &= \mqty*(1+2\beta & \alpha - \beta & 0 & \cdot & \cdot & 0 \\ -\alpha -\beta & 1+2\beta & \alpha-\beta & 0 & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & \cdot & 0 & -\alpha -\beta & 1+2\beta & \alpha-\beta \\ 0 & \cdot & \cdot & 0 & -\alpha-\beta & 1+2\beta);\\
	R &= \mqty*(1-2\beta & -\alpha + \beta & 0 & \cdot & \cdot & 0 \\ \alpha +\beta & 1-2\beta & -\alpha+\beta & 0 & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & \cdot & 0 & \alpha +\beta & 1-2\beta & -\alpha+\beta \\ 0 & \cdot & \cdot & 0 & \alpha+\beta & 1-2\beta);\\
	\vb{b} &= \mqty*((\alpha + \beta)(u_0^n + u_0^{n+1}) \\ 0 \\ \cdot \\ 0 \\ (-\alpha + \beta)(u_M^n + u_M^{n+1})) \equiv \frac{r}{2}\mqty*((\alpha + \beta)(g_0(t_n) + g_0(t_{n+1})) \\ 0 \\ \cdot \\ 0 \\ (-\alpha + \beta)(g_M(t_n) + g_M(t_{n+1}))).
\end{align}
Notiamo poi subito che le matrici $L$ ed $R$ possono essere rispettivamente scritte in forma $(I-A)$ e $(I+A)$ dove $A$ è definita come
\begin{equation}
	A = \mqty*(-2\beta & -\alpha + \beta & 0 & \cdot & \cdot & 0 \\ \alpha +\beta & -2\beta & -\alpha+\beta & 0 & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & \cdot & 0 & \alpha +\beta & -2\beta & -\alpha+\beta \\ 0 & \cdot & \cdot & 0 & \alpha+\beta & -2\beta). 
\end{equation}
Ottenendo quindi
\begin{equation}
	\left(I-A\right)\vb{u}^{n+1} = \left(I+A\right)\vb{u}^n + \vb{b}
\end{equation}

\subsubsection{Note extra riguardanti Crank-Nicolson}
Come abbiamo visto nelle sezioni precedenti, molto spesso si ha a che fare con una equazione differenziale nella forma
\begin{equation}
	u_t = {(\alpha(x,t)u_x)}_x + \beta(x,y)u_x + \gamma(x,y)u
\end{equation}
è possibile in questo caso compiere la seguente discretizzazione per la costruzione dell'integratore di Crank-Nicolson:
\begin{align}
	u_t &\Rightarrow \frac{1}{k}\delta_t U_m^n,\\
	{(\alpha(x,t)u_x)}_x &\Rightarrow \frac{1}{2h}\left(\alpha_{m+\frac{1}{2}}^n \frac{\delta_x U_m^n}{h} - \alpha_{m-\frac{1}{2}}^n \frac{\delta_x U_{m-1}^n}{h}\right) + \frac{1}{2h}\left(\alpha_{m+\frac{1}{2}}^{n+1} \frac{\delta_x U_m^{n+1}}{h} - \alpha_{m-\frac{1}{2}}^{n+1} \frac{\delta_x U_{m-1}^{n+1}}{h}\right)\\
	\beta(x,t)u_x &\Rightarrow \frac{1}{4h}(\beta_m^n(U_{m+1}^n - U_{m-1}^n) + \beta_m^{n+1}(U_{m+1}^{n+1} - U_{m-1}^{n+1}))\\
	\gamma(x,t)u &\Rightarrow \frac{1}{2}(\gamma_m^n U_m^n + \gamma_m^{n+1} U_m^{n+1})
\end{align}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Oscillatore armonico smorzato con rumore stocastico}
Consideriamo un oscillatore armonico immerso in un fluido soggetto a due tipi di forze, una stocastica dovuta a urti casuali con le molecole del fluido ed una sistematica dovuta alla resistenza che il fluido oppone al suo moto. Abbiamo l'equazione di Langevin:
\begin{align}
	\dot{q} &= p\\
	\dot{p} &= -\omega^2 q - \gamma p + \sqrt{2T\gamma}\,\,\xi(t)
\end{align}
Dalle considerazioni fatte nella sezione \ref{sec:fluttuazione_dissipazione}, possiamo calcolare la distribuzione di equilibrio
\begin{equation}
	\rho(q,p) = Z^{-1} \exp{-\frac{\frac{p^2}{2}+\frac{\omega^2q^2}{2}}{T}}
\end{equation}
dove in questo caso \(Z^{-1}=\omega/2T\pi\).

Se vogliamo invece trattare il problema utilizzando le variabili azione-angolo standard per un oscillatore armonico ideale
\begin{equation}
	J = \frac{p^2}{2\omega} + \frac{\omega q^2}{2} \quad\quad\quad\quad \theta = \atan\left(\frac{\omega q}{p}\right)
\end{equation}
ci troviamo ad avere
\begin{equation}
	\rho(J, \theta) = \frac{\omega}{2\pi T} \exp\left(-\frac{J\omega}{T}\right)
\end{equation}
dove figura immediatamente come \(\theta\) sia uniformemente distribuita all'equilibrio (e che quindi compare come un termine \(1/2\pi\) nella costante di normalizzazione).\\
Portando avanti i conti e sfruttando in particolare
\begin{equation}
	q = \sqrt{\frac{2J}{\omega}}\,\,\sin{\theta} \quad\quad\quad p = \sqrt{2J\omega}\,\,\cos{\theta}
\end{equation}
si ha che
\begin{align}
	\dot{J} &= \frac{p\dot{p}}{\omega} + \omega q\dot{q}\\
	&= -2J\gamma \cos^2\theta \,+\, \xi(t)\sqrt{\frac{4JT\gamma}{\omega}}\cos\theta \\
	\dot{\theta} &= \frac{\omega p \dot{q} - \omega q \dot{p}}{\omega^2 q^2 + p^2}\\
	&= \frac{\omega p^2 -\omega q(-\omega^2 q - \gamma p + \sqrt{2 T \gamma}\xi(t))}{\omega^2 q^2 + p^2}\\
	&= \omega + \gamma \sin{\theta}\cos{\theta} - \sqrt{\frac{T\gamma}{J\omega}}\xi(t)\sin{\theta}
\end{align}
per ottenere una equazione descrivente l'evoluzione di \(\rho\), è possibile ragionare in termini di approssimazione di media sulla variabile angolare: se imponiamo che
\begin{equation}
	\frac{\gamma}{\omega} \ll 1
\end{equation}
possiamo considerare i tempi di evoluzione della variabile angolare \(\theta \) sufficientemente piccoli rispetto ai tempi di evoluzione della variabile azione \(J\) e trattare quindi il sistema con il metodo dell'averaging principle trattato in sezione~\ref{subsec:averaging_principle}. L'applicazione del principio equivale a trattare \(\theta \) come distribuita uniformemente lungo la superficie toroidale, pur essendo comunque i suoi tempi effettivi di rilassamento a tale distribuzione comparabili a quelli della variabile azione. 

Scrivendo in prima istanza l'equazione di evoluzione di $\rho$:
\begin{multline}
	\pdv{\rho}{t} = \gamma \pdv{}{J} 2J\cos^2\theta\,\,\rho + \pdv{}{J} J\frac{2T\gamma}{\omega}\cos^2\theta \pdv{}{J} \rho \\+ \pdv{}{\theta} \left(\omega + \gamma \sin{\theta}\cos{\theta}\right) \rho + \pdv{}{\theta} \frac{T\gamma}{J\omega}\sin^2\theta \pdv{}{\theta}\rho
\end{multline} 
possiamo quindi trascurare le derivate angolari e considerare $\theta$ come uniformemente distribuita lungo la superficie toroidale, ottenendo quindi
\begin{equation}
	\pdv{\rho_\theta}{t} = \gamma \pdv{}{J} J\rho_\theta + \pdv{}{J} J\frac{T\gamma}{\omega} \pdv{}{J} \rho_\theta
	\label{eq:pendolo_fokker_planck_1}
\end{equation}
dove abbiamo sfruttato \(\frac{1}{2\pi}\int_0^{2\pi}\cos^2\theta\, d\theta = \frac{1}{2}\).\\

\subsection{Analisi numerica}

Si cerca ora un riscontro della teoria nei risultati numerici. Dove non indicato diversamente, vengono assegnati i valori numerici riportati in Tabella~\ref{tab:valori_1}.

\begin{table}[!h]
	\centering
	\begin{tabular}{cS}
	\toprule
	Costante & {Valore assegnato} \\
	\midrule
	$\omega$	& 1.	\\
	$\gamma$	& 0.01	\\
	$T$			& 10.	\\
	$dt$		& 0.1	\\
	\bottomrule
	\end{tabular}
	\caption{Valori numerici usati per le costanti del modello.}
	\label{tab:valori_1}
\end{table}
Come condizioni iniziali, viene considerata per la variabile azione $I$ una distribuzione di probabilità iniziale gaussiana a media pari a 10 e varianza unitaria. Per la variabile angolare $\theta$ vengono trattati separatamente il caso in cui la distribuzione iniziale è uniforme (quindi già rilassata) ed il caso in cui invece la distribuzione iniziale è $\rho_0(\theta)=\delta(\theta_0)$.

\subsubsection{Rilassamento all'equilibrio}

Dati i valori in Tabella~\ref{tab:valori_1}, abbiamo il grafico della distribuzione di equilibrio per \(\rho(J)\) in Fig.~\ref{fig:equilibrio_J_1}, mentre per \(\rho(\theta)\) si ha una distribuzione uniforme lungo l'intervallo \([0,2\pi]\).

\begin{figure}[p]
	\centering
	\includegraphics[width=\textwidth]{selected/equilibrium.jpg}
	\caption{Distribuzione di equilibrio per \(\rho(I)\), la variabile angolare ed il suo parametro di normalizzazione $1/2\pi$ non vengono considerati.}
	\label{fig:equilibrio_J_1}
\end{figure}

Possiamo in primo luogo simulare numericamente l'equazione di Langevin usando l'integratore descritto in Sezione~\ref{sec:integratore} ed ottenere i grafici riportati in Figura~\ref{fig:langevin_0} e~\ref{fig:langevin_150} per quel che riguarda la condizione inziale con variabile angolo concentrata su di un unico valore ed i grafici riportati in Figura~\ref{fig:langevin_uniform_0} e~\ref{fig:langevin_uniform_150} per quel che riguarda la condizione inziale con variabile angolo uniformemente distribuita. Oltre a questi grafici sono anche stati prodotti dei video che permettono di apprezzare l'evoluzione delle distribuzioni. Per ogni grafico e video prodotto sono stati considerati in blocco 100000 diverse realizzazioni del sistema perturbato.

Le diverse condizioni iniziali sulla variabile angolare non interferiscono con il processo di rilassamento della variabile azione. Si osserva inoltre come la variabile angolare risulta avere tempi di rilassamento superiori a quelli della variabile azione.

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/langevine_t_0.jpg}
	\caption{Integrazione dell'equazione di Langevin, condizioni iniziali con variabile angolo piccata su di un unico valore. La curva arancione rappresenta la distribuzione di equilibrio.}
	\label{fig:langevin_0}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/langevine_t_150.jpg}
	\caption{Integrazione dell'equazione di Langevin, rilassamento della variabile azione, condizioni iniziali con variabile angolo piccata su di un unico valore. La curva arancione rappresenta la distribuzione di equilibrio.}
	\label{fig:langevin_150}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/langevine_uniform_t_0.jpg}
	\caption{Integrazione dell'equazione di Langevin, condizioni iniziali con variabile già rilassata a distribuzione uniforme. La curva arancione rappresenta la distribuzione di equilibrio.}
	\label{fig:langevin_uniform_0}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/langevine_uniform_t_150.jpg}
	\caption{Integrazione dell'equazione di Langevin, rilassamento della variabile azione, condizioni iniziali con variabile già rilassata a distribuzione uniforme. La curva arancione rappresenta la distribuzione di equilibrio.}
	\label{fig:langevin_uniform_150}
\end{figure}
Per quanto riguarda invece l'equazione differenziale \eqref{eq:pendolo_fokker_planck_1}, vogliamo verificare che, con equivalenti condizioni iniziali (distribuzione gaussiana a media 10 e varianza unitaria), questa descrivi correttamente il rilassamento della variabile azione fino alla distribuzione di equilbrio, descritto dalla integrazione numerica dell'equazione di Langevin. Per integrarla numericamente si fa uso dello schema implicito di Crank-Nicolson impostando come parametri ulteriori necessari\footnote{\(L\) rappresenta la regione considerata per l'integrazione, si pone una barriera assorbente per valori di \(J\) superiori al limite considerato, questo inevitabilmente porterà su lunghi periodi a perdita di flusso, tuttavia, ponendo la barriera molto in avanti, dove la distribuzione teorica di equilibrio assume valori prossimi allo zero, si può ritenere questa perdita trascurabile.}:
\begin{table}[!]
	\centering
	\begin{tabular}{cc}
	\toprule
	Costante & {Valore assegnato} \\
	\midrule
	\(L\)		& \([0,200]\)	\\
	\(dx\)		& \(0.05\)	\\
	\bottomrule
	\end{tabular}
	\caption{Parametri extra per integratore di Crank-Nicolson.}
	\label{tab:valori_2}
\end{table}
L'integrazione di Crank, sovrapposta agli istogrammi ottenuti dall'integrazione dell'equazione di Langevin e alla curva della distribuzione stazionaria, è rappresentata nelle Figure~\ref{fig:crank_1},~\ref{fig:crank_2} e~\ref{fig:crank_3} sovrapponendo i dati di Langevin con distribuzione angolare piccata su di un unico valore e nelle Figure~\ref{fig:crank_4},~\ref{fig:crank_5} e~\ref{fig:crank_6} con i dati di Langevin con distribuzione angolare uniforme. Si osserva in ogni momento perfetta coerenza tra i risultati dei due integratori.

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_t_0.jpg}
	\caption{Integratore di Crank-Nicolson, condizioni iniziali per i dati su istogramma con variabile angolo piccata su di un unico valore.}
	\label{fig:crank_1}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_t_25.jpg}
	\caption{Integratore di Crank-Nicolson, tempo intermedio, condizioni iniziali per i dati su istogramma con variabile angolo piccata su di un unico valore.}
	\label{fig:crank_2}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_t_150.jpg}
	\caption{Integratore di Crank-Nicolson, rilassamento, condizioni iniziali per i dati su istogramma con variabile angolo piccata su di un unico valore.}
	\label{fig:crank_3}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_uniform_t_0.jpg}
	\caption{Integratore di Crank-Nicolson, condizioni iniziali per i dati su istogramma con variabile angolo uniformemente distribuita.}
	\label{fig:crank_4}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_uniform_t_25.jpg}
	\caption{Integratore di Crank-Nicolson, tempo intermedio, condizioni iniziali per i dati su istogramma con variabile angolo uniformemente distribuita.}
	\label{fig:crank_5}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=0.8\textwidth]{selected/crank_uniform_t_150.jpg}
	\caption{Integratore di Crank-Nicolson, rilassamento, condizioni iniziali per i dati su istogramma con variabile angolo uniformemente distribuita.}
	\label{fig:crank_6}
\end{figure}


\end{document}