Un algoritmo di Monte Carlo per la simulazione di un …...un algoritmo di Monte Carlo e riduce...

UNIVERSITÀ DEGLI STUDI DI PISA

DIPARTIMENTO DI FISICA “E. FERMI”

CORSO DI LAUREA IN FISICA

UN ALGORITMO DI MONTE CARLO PER LASIMULAZIONE DI UN SISTEMA FRUSTRATO

TESI DI LAUREA TRIENNALE

Presentata da:GUIDO MASELLA

Relatore:Dott.

GIANCARLO CELLA

Controrelatore:Prof.

ROSA POGGIANI

Anno Accademico 2012 – 2013

2

Indice

1 Teoria Generale 31.1 Metodi di Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.1.1 Critical slowing-down . . . . . . . . . . . . . . . . . . . . . . . . 41.2 Modello di Ising bidimensionale . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.1 Heat Bath e algoritmi a mosse locali. . . . . . . . . . . . . . . . . 61.2.2 Swendsen Wang . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Generalizzazioni del modello di Ising e modello O(N) . . . . . . . . . . . 8

2 L’algoritmo 112.1 Azione Standard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.1.1 L’update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2 Azione Symanzik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2.1 Studio di un caso particolare: la regola H2 . . . . . . . . . . . . . . 192.2.2 Studio delle regole: classi ed applicabilità . . . . . . . . . . . . . . 20

2.2.2.1 Classi di regole e altre semplificazioni . . . . . . . . . . 202.3 Implementazione dell’algoritmo . . . . . . . . . . . . . . . . . . . . . . . 25

3 Risultati 293.1 Misure da eseguire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2 Analisi dell’applicabilità delle regole. . . . . . . . . . . . . . . . . . . . . 363.3 Energia: confronto con l’Heat Bath. . . . . . . . . . . . . . . . . . . . . . 363.4 Critical Slowing-Down . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

i

ii INDICE

Elenco delle figure

1.2.1 Esempio update Swendsen e Wang. . . . . . . . . . . . . . . . . . . . . . 9

2.1.1 Rappresentazione fattore di Boltzmann. . . . . . . . . . . . . . . . . . . . 122.1.2 Espansione del fatttore di Boltzmann. . . . . . . . . . . . . . . . . . . . . 122.1.3 Quattro passi successivi algoritmo SW standard. . . . . . . . . . . . . . . . 142.2.1 Placchetta di tre siti per l’azione Symanzik. . . . . . . . . . . . . . . . . . 152.2.2 Fattore di Boltzmann hamiltoniana Symanzik. . . . . . . . . . . . . . . . . 162.2.3 Espansione fattore di Boltzmann Symanzik . . . . . . . . . . . . . . . . . 172.2.4 Regole di update algoritmo Symanzik. . . . . . . . . . . . . . . . . . . . . 182.3.1 Diagramma delle classi che implementano l’algoritmo di update. . . . . . . 272.3.2 classiDiagramma delle classi utili all’acquisizione di statistica. . . . . . . . 28

3.2.1 Applicabilità regole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.2 Applicabilità regole, classe 1. . . . . . . . . . . . . . . . . . . . . . . . . . 323.2.3 Applicabilità regole, classe 2. . . . . . . . . . . . . . . . . . . . . . . . . . 333.2.4 Applicabilità regole, classe 3. . . . . . . . . . . . . . . . . . . . . . . . . . 343.2.5 Applicabilità regole, classe 4. . . . . . . . . . . . . . . . . . . . . . . . . . 353.3.1 Energia regola H1, a = 1/2. . . . . . . . . . . . . . . . . . . . . . . . . . 373.3.2 Energia regola H2, a = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.3.3 Energia regola D2, a = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.4.1 FSS ansatz, regola H2 a = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 403.4.2 FSS ansatz, regola H1 a = 1/2. . . . . . . . . . . . . . . . . . . . . . . . . 413.4.3 FSS ansatz, regola D2 a = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 42

iii

iv ELENCO DELLE FIGURE

Introduzione

Per valutare il valore di aspettazione di osservabili su sistemi fisici a molti gradi di libertàdescritti con i metodi della meccanica statistica, bisogna ricorrere a metodi numerici e inparticolare a metodi di Monte Carlo. In generale la struttura di questi metodi è semprela stessa: data la distribuzione di equilibrio per le configurazioni del sistema si genera daquesta un insieme molto esteso di sistemi campione su cui verranno determinati i valori diaspettazione delle osservabili cercate.

Per sistemi molto complessi, con molti gradi di libertà si procede costruendo un opportu-no processo stocastico avente come distribuzione di equilibrio la distribuzione di equilibriodel sistema fisico da simulare. In questo modo le medie temporali sul processo stocasticoconvergeranno alle medie sulla distribuzione di equilibrio. In generale tale processo stoca-stico sarà una catena di Markov la cui matrice di transizione deve avere un unico autovaloreunitario, corrispondente alla distribuzione di equilibrio. Ciò si può ottenere semplicementeimponendo il bilancio dettagliato e la ergodicità.

Un problema di questi metodi risiede nel fatto che stati successivi della catena di Markovpossono essere fortemente correlati riducendo così l’efficacia dell’importance sampling, ela riduzione si può quantificare in termini di un tempo caratteristico di correlazione. Inalcuni regimi tale parametro può divergere, in particolare in prossimità di un punto critico(il cosiddetto critical slowing down), oppure a bassa temperatura in presenza di frustrazione.

Per quanto riguarda la costruzione di algoritmi per modelli di spin possiamo distingueredue classi fondamentali. La prima classe comprende algoritmi con metodi di update localicome quelli proposti da Metropolis o algoritmi come l’Heat Bath. Caratteristica comunedegli algoritmi di questa prima classe è che tutti i casi presenti in letteratura presentano unsevero slowing down. Gli algoritmi della seconda classe, che prevede metodi di update col-lettivi (non locali), possono ridurre drasticamente il critical slowing down. Da questo pun-to di vista estremamente efficienti sono gli algoritmi basati su un metodo originariamenteproposto da Swendsen e Wang.

In questa tesi ci focalizzeremo in particolare sull’applicazione di tali algoritmi ad unaclasse di modelli detti modelli sigmaO(n) non lineari, nel limite del continuo. Se il modelloè descritto da una azione non minimale il sistema statistico a cui viene applicato l’algoritmodi Swendsen Wang risulta parzialmente frustrato, e l’efficienza tende a peggiorare.

Verrà dapprima discusso ed implementato l’algoritmo Heat Bath che sarà utilizzato

1

2 ELENCO DELLE FIGURE

come termine di confronto per gli altri algoritmi analizzati.Verrà quindi proposta e discussa una applicazione del metodo di Swendsen e Wang ad

un sistema descritto da una Hamiltoniana standard (modello di Ising).Successivamente verrà discussa la costruzione di possibili varianti dell’algoritmo alla

Swendsen e Wang utilizzato in [3] per sistemi descritti da una Hamiltoniana modificatasecondo un metodo suggerito da K. Symanzik.

Le proprietà di tali soluzioni verranno studiate numericamente, cercando di determinarei fattori che ne influenzano il comportamento critico, anche alla luce delle informazionipresenti in letteratura.

Capitolo 1

Teoria Generale

1.1 Metodi di Monte CarloIl problema da risolvere è quindi quello di generare, data una certa distribuzione (π) diequilibrio del sistema considerato, un insieme molto esteso di configurazioni campione sucui verranno determinati i valori di aspettazione delle osservabili cercate. Principalmentei metodi di Monte Carlo si differenziano in metodi statici e metodi dinamici. I metodidi Monte Carlo statici consistono nel creare una sequenza statisticamente indipendente disistemi campione a partire dalla distribuzione data. Per sistemi molto complessi, con unnumero elevato di gradi di libertà, dove i metodi statici risultano praticamente infattibili,si procede costruendo un opportuno processo stocastico avente π come distribuzione diequilibrio. In generale tale processo stocastico sarà una catena di Markov. Bisogna oradeterminare le le proprietà generali della catena di Markov.

In generale definiamo una catena di Markov su uno spazio discreto di stati S come unasuccessione di stati X0, X1, X2, · · · ∈ S tale che transizioni successive Xt → Xt+1 sia-no statisticamente indipendenti, ovvero tale che lo stato Xt+1dipenda esclusivamente dallostato Xt e non dagli stati precedenti. Bisogna quindi assicurarsi che il processo raggiungal’equilibrio, che la sua distribuzione di equilibrio sia proprio la distribuzione π data e che lemedie temporali sul processo stocastico convergano alle medie sulla distribuzione di equi-librio (ergodicità). Sia px→y = P(Xt+1 = y |Xt = x) la probabilità di transizione dallostato x ∈ S allo stato y ∈ S definiamo quindi P = {px→y}x,y∈S la matrice di transizionedella catena di Markov e p(n)x→y = P(Xt+n = y |Xt = x) come la probabilità di transizioneda x a y in n passi. Ai fini del nostro problema richiediamo quindi che la catena di Markovsoddisfi le seguenti condizioni:

• Ergodicità: Per ogni x, y ∈ S esiste n > 0 tale che p(n)x→y > 0

• Stazionarietà:∀y ∈ S

∑x

πxpx→y = πy (1.1.1)

3

4 CAPITOLO 1. TEORIA GENERALE

Nella maggior parte dei casi conviene imporre una condizione più forte della stazionarietàma più semplice da controllare:

• Bilancio dettagliato: Per ogni x, y ∈ S,

πxpx→y = πypx→y (1.1.2)

Si può dimostrare (vedere [8]) che queste condizioni assicurano che qualsiasi sia lo stato dipartenza, per t → ∞ il processo converge all’equilibrio e le medie temporali convergonoalle medie sulla distribuzione π di equilibrio.

1.1.1 Critical slowing-downIl problema dei metodi dinamici di Monte Carlo risiede nel fatto che stati successivi del pro-cesso stocastico considerato ovvero, nel nostro caso della catena di Markov, sono correlati,magari anche fortemente.

Si vuole ora capire quanto e come la correlazione tra stati successivi della catena di Mar-kov influisca sulla misura di una certa osservabile f definita sul sistema da noi considerato.Sia quindi {ft} un processo stocastico stazionario a valori reali con media

µ = 〈ft〉 (1.1.3)

definiamo la funzione di autocorrelazione non normalizzata come

Cff (t) ≡ 〈fsfs+t〉 − µ2 (1.1.4)

e ρ(t) la funzione di autocorrelazione normalizzata:

ρff (t) ≡C(t)

C(0). (1.1.5)

Tipicamente ρ(t) decade come e−|t|/τ per grandi t. Possiamo definire due parametri carat-teristici: un tempo di autocorrelazione esponenziale

τexp,f = lim supt→∞

t

− ln |ρff (t)|, (1.1.6)

che rappresenta il tempo di rilassamento del modo più lento nel sistema, ed un tempo diautocorrelazione integrato

τint,f =1

2

∞∑t=−∞

ρff (t) =1

2+∞∑t=1

ρff (t) (1.1.7)

1.2. MODELLO DI ISING BIDIMENSIONALE 5

che invece controlla l’errore statistico nella misura di Monte Carlo di 〈f〉:

V ar[f̄ ] =1

n2

n∑r,s=1

Cff (r − s) =1

n

n−1∑t=−(n−1)

Cff (t) ≈1

n2τint,fCff (0) (1.1.8)

dove f̄ è la media

f̄ =1

n

n∑t=1

ft (1.1.9)

Ne segue che la varianza della media è più grande di un fattore 2τint,f di quanto sarebbestata se gli ft fossero stati statisticamente indipendenti tra loro e che quindi il numero dicampioni “effettivamente indipendenti” in un run di lunghezza n sono n

2τint,f. Ai fini della

stima dei valori di aspettazione dell’osservabile f possiamo dire che τexp,f rappresenta l’in-formazione sul transiente iniziale di misure da scartare mentre il tempo di autocorrelazioneintegrato τint,f rappresenta l’informazione sulla lunghezza effettiva del run. In alcuni regi-mi, in prossimità di un punto critico o a basse temperature, in presenza di frustrazione, iltempo di autocorrelazione τint,f può divergere, in particolare vicino ad un punto critico iltempo di autocorrelazione diverge come

τ ∼ ξz (1.1.10)

dove ξ è la lunghezza di correlazione del sistema considerato e z è un esponente criticodinamico. Questo fenomeno, noto in letteratura come critical slowing-down riduce note-volmente il numero di campioni statisticamente indipendenti in un run di lunghezza data diun algoritmo di Monte Carlo e riduce così l’efficacia dell’importance sampling.

1.2 Modello di Ising bidimensionaleConsideriamo inizialmente il modello di Ising per un reticolo bidimensionale quadrato inassenza di campo esterno. Se il reticolo è finito imponiamo condizioni periodiche al con-torno in modo da poter trascurare effetti di superficie. Per ogni sito del reticolo indichiamocon n la posizione del sito, con il vettore unitario µ ∈ {ex, ey} una delle due direzionipossibili sul reticolo e con σn una variabile discreta definita sul sito che può prendere comevalori solo ±1. Consideriamo quindi l’hamiltoniana del modello di Ising:

HIsingσ = −

∑n,µ

Jn,n+µσnσn+µ (1.2.1)

dove con n + µ si indica il sito primo vicino al sito n lungo la direzione µ e Jn,n+µ ≥0 rappresenta l’accoppiamento ferromagnetico tra i siti n e n + µ. Per un modello di


Ising generalmente Jn,n+µ è una costante indipendente dalla posizione sul reticolo. Dallameccanica statistica sappiamo che la distribuzione di equilibrio π del reticolo è data da

π(σ) =1

ZBIsingσ =

1

Ze−βH

Isingσ (1.2.2)

dove Z è la funzione di partizione del sistema.

1.2.1 Heat Bath e algoritmi a mosse locali.Gli algoritmi volti a simulare la distribuzione di equilibrio (1.2.2) si possono distinguerein due classi principali; la prima classe comprende algoritmi con metodi di update sostan-zialmente locali, ovvero fissato un sito n del reticolo si propone di cambiare il valore dellavariabile associata σn. Spesso si può fare in modo che il “nuovo” valore della variabile σnsia indipendente dal valore precedente. Un possibile algoritmo di questa classe è l’HeatBath, dove la probabilità condizionata di σn, dati tutti gli altri spin {σn′} con n′ 6= n è

Pr[σn

∣∣∣{σn}n′ 6=n

]= const

({σn}n′ 6=n

)× exp

(−βσn

∑η

ση

)(1.2.3)

e la somma è su tutti i siti η primi vicini di n. L’algoritmo prevede inoltre che tale opera-zione sia ripetuta su tutti i siti del reticolo in maniera casuale o meno1. Nel caso dell’HeatBath il nuovo valore di σn è effettivamente indipendente dal valore precedente e dipende,a parte una costante, dallo spin dei primi vicini ovvero dall’interazione che i primi viciniesercitano sul sito stesso.

Questa classe di algoritmi è affetta da un severo slowing down, in particolare vicinoad una transizione di fase del secondo ordine il tempo di autocorrelazione diverge, comeabbiamo già visto, con una legge del tipo

τ ∼ min(L, ξ)z (1.2.4)

dove L è la dimensione lineare del sistema e ξ è la lunghezza di correlazione del sistemae z è un esponente critico che per questi algoritmi assume valori tipicamente vicini a 2.Questo è dovuto principalmente al fatto che l’informazione sulla nuova configurazione ètrasmessa da un sito solo ai suoi primi vicini, per ogni step dell’algoritmo, e affinché ilsistema evolva in una situazione completamente nuova e scorrelata dalla precedente, biso-gna che quest’informazione raggiunga distanze dell’ordine di ξ ovvero della lunghezza dicorrelazione. Siccome possiamo immaginare che l’informazione propagandosi esegua unrandom walk nel reticolo, ci aspettiamo che il tempo necessario per percorrere una distanzaξ sia τ ∼ ξ2 e questo in qualche modo giustifica il valore presente in letteratura [10, 11] peralgoritmi a mosse locali come l’Heat Bath.

1Ad esempio si può partizionare i siti come una scacchiera, e alterare sequenzialmente le caselle “bianche”e quelle “nere”.

1.2. MODELLO DI ISING BIDIMENSIONALE 7

1.2.2 Swendsen WangUna seconda classe di algoritmi viene introdotta proprio per risolvere il severo slowing-down intrinseco degli algoritmi a mosse locali. Gli algoritmi di questa classe sono quin-di algoritmi con metodi di update collettivi (non locali) e ne esistono di diversi tipi. Inquesta tesi tratteremo solo algoritmi basati su un idea originale sviluppata da Swendsen eWang [15]. La soluzione proposta consiste nel “complicare” il modello aggiungendo nuovevariabili e quindi simulare questo modello aumentato.

Consideriamo la distribuzione di equilibrio del modello di Ising (1.2.2). Notiamo che lafunzione di partizione è, in questo caso, solo una costante di normalizzazione e possiamotrascurarla. Interessiamoci solo al fattore di Boltzmann BIsing

σ relativo ad una configurazioneσ dell’intero reticolo e scriviamolo come prodotto

BIsingσ =

∏n,µ

BIsingn,µ (1.2.5)

dove ciascun termine sarà della forma

BIsingn,µ = eβJn,n+µδσn,σn+µ + e−βJn,n+µ δ̄σn,σn+µ (1.2.6)

dove δA,B = 1 se A = B altrimenti δA,B = 0 e δ̄A,B = 1 − δA,B. Per ogni fattore Bn,µintroduciamo una nuova variabile ln,µ ∈ {0, 1} e scriviamo

B̂Isingn,µ =

(δσn,σn+µ , δ̄σn,σn+µ

)·(W11 W12

W21 W22

)·(ln,µl̄n,µ

)(1.2.7)

dove abbiamo definito l̄n,µ = 1− ln,µ. Le nuove variabili ln,µ appena introdotte corrispon-dono a dei collegamenti tra i siti n e n + µ: se ln,µ = 1 allora i due siti sono collegati e laloro orientazione relativa non può cambiare; al contrario, come precedentemente imposto,se ln,µ = 0 allora gli spin σn e σn+µ sono completamente disaccoppiati tra loro.

Il metodo di update di un algoritmo alla Swendsen e Wang consiste quindi nel sud-dividere il reticolo in cluster di siti collegati tra loro e quindi “ribaltare” collettivamenteciascun cluster indipendentemente con probabilità 1

2.

Per specificare completamente l’algoritmo bisogna ora determinare i coefficienti Wij

presenti nella (1.2.7). Per consistenza bisogna imporre che, sommando su tutte le possi-bili configurazioni di ln,µ si ottenga nuovamente la distribuzione (1.2.2), tale condizioneequivale ad imporre che ∑

ln,µ

B̂Isingn,µ = kn,µ(σ)Bn,µ (1.2.8)

da cui si ottiene che

(W11 +W12) e−βJn,n+µ = (W22 +W21) e

+βJn,n+µ , (1.2.9)


inoltre richiediamo che quando ln,µ = 0 gli spin σn e σn+µ siano disaccoppiati, da cui

W22 = W12 (1.2.10)

e che quando ln,µ = 1 le relative orientazioni rimangano “congelate”, ovvero

W11W21 = 0. (1.2.11)

Abbiamo quindi due modi di scrivere il fattore di Boltzmann del modello aumentato:

B̂Isingn,µ =


)·(e2βJn,n+µ − 1 1

0 1

)·(ln,µl̄n,µ

)(1.2.12)

B̂Isingn,µ =


)·(

0 11− e2βJn,n+µ 1

)·(ln,µl̄n,µ

)(1.2.13)

ma tra queste, per un accoppiamento ferromagnetico, la seconda non è accettabile in quantotutti gli gli elementi della matrice in (1.2.13) non sono dello stesso segno, e quindi nonpossono essere interpretati come probabilità a parte un fattore di normalizzazione comune.La distribuzione di probabilità dei ln,µ condizionata dalla configurazione degli spin è quindidata dalle relazioni

Pr[`n,µ = 1

∣∣δσn,σn+µ = 1]

=W11

W11 +W12

= 1− e−2βJn,n+µ (1.2.14)

Pr[`n,µ = 0


=W12

W11 +W12

= e−2βJn,n+µ (1.2.15)

Pr[`n,µ = 1


=W21

W21 +W22

= 0 (1.2.16)

Pr[`n,µ = 0


=W22

W21 +W22

= 1 (1.2.17)

1.3 Generalizzazioni del modello di Ising e modello O(N)

Al posto di definire, per ogni sito del reticolo, una variabile discreta σn che può assumeresolo i valori ±1 come nel modello di Ising, possiamo definire un vettore a due componenti

1.3. GENERALIZZAZIONI DEL MODELLO DI ISING E MODELLO O(N) 9

Figure 1.2.1: Semplice schema che mostra il procedimento di update di un algoritmo allaSwendsen e Wang per un reticolo 6 × 6. I punti blu e rossi rappresentano i valori possibilidella variabile σn mentre i tratti rappresentano i link tra i siti. Lo schema consiste quindinella costruzione del cluster e nel suo successivo ribaltamento. La probabilità che si formiun link tra due siti con lo stesso orientamento è data dall’Equazione (1.2.14).

pn che può assumere un certo numero intero q di valori uniformemente distribuiti su unacirconferenza:

θn =2πh

q, h = 1, . . . , q (1.3.1)

dovepn ≡ (cos θn, sin θn). (1.3.2)

L’hamiltoniana di un tale modello sarà quindi

HPotts = −∑n,µ

Jn,n+µ(pn · pn+µ

)(1.3.3)

e si riduce al modello di Ising per q = 2.Una ulteriore generalizzazione consiste nel considerare il limite del modello di Potts per

q → ∞ questo è il caso del modello planare XY [6, 9] dove ad ogni sito è associata unavariable θ̄n continua a valori nell’intervallo [0, 2π]. L’hamiltoniana del modello è

HXY = −∑n,µ

Jn,n+µ (rn · rn+µ) (1.3.4)

dove rn è un vettore sul piano di lunghezza unitaria, ovvero

rn = (cos θ̄n, sin θ̄n) (1.3.5)

Più in generale si possono considerare le variabili di spin non come vettori in uno spaziobidimensionale ma come vettori di lunghezza unitaria in uno spazio di dimensione N qual-siasi; questo è il caso dei cosiddetti modelli O(N) [13]. Un primo esempio è dato dallahamiltoniana standard

HSTD = −∑n,µ

φn · φn+µ (1.3.6)


dove i φn sono vettori a N componenti normalizzati in modo da avere φn · φn = 1.In questa tesi considereremo un altro modello sempre appartenente alla classe di modelli

O(N). L’hamiltoniana è data in questo caso da

HSYM = −∑n,µ

(4

3φn · φn+µ −

1

12φn · φn+2µ

)(1.3.7)

è stata proposta inizialmente da K. Symanzik e da lui prende il nome. Accenniamo breve-mente alle motivazioni che stanno dietro alla (1.3.7). L’idea di base è che si vuole studiareuna certa teoria di campo bidimensionale, e per farlo si approssima lo spazio discretizzan-dolo su un reticolo. Gli operatori differenziali presenti nella hamiltoniana non discretizzatavengono approssimati con differenze finite, e nel fare questo si introducono errori detti ar-tefatti di reticolo. Sia la (1.3.6) che la (1.3.7) possono essere viste come approssimazionidella stessa teoria di campo nel continuo. La differenza è che, detto a il passo reticolare,l’errore dovuto alla discretizzazione è O(a) nel primo caso, O(a2) nel secondo.

Notare che, al contrario della (1.3.6), la (1.3.7) include anche un termine di interazionecon terzi vicini (di natura antiferromagnetica).

Capitolo 2

L’algoritmo

2.1 Azione StandardConsideriamo ora l’azione standard di un modello O(N). Come già abbiamo detto l’hamil-toniana è della forma

HSTD = −∑n,µ

φn · φn+µ (2.1.1)

con φn·φn = 1. L’applicazione di un algoritmo alla Swendsen e Wang al modello sigma nonlineare necessita di una procedura di embedding in grado di ridurre il sistema statistico adun modello di Ising con interazioni casuali. Lo schema di embedding proposto da Wolff [5]consiste nell’introdurre, per ogni sito, una variabile di Ising σn ∈ {1,−1} tale che, sceltoun vettore unitario r nello spazio N -dimensionale dei φn valga la decomposizione

φn = φ⊥n + σn |φn · r|r (2.1.2)

Il vettore r è scelto casualmente. In questo modo il ribaltamento di uno spin (σn → −σn)è equivalente ad una riflessione di φn nel piano perpendicolare alla direzione r. La partedell’hamiltoniana (2.1.1) dipendente da σ = {σn} diventa quindi

HSTDσ = −

∑n,µ

Jn,n+µσn · σn+µ (2.1.3)

del tutto equivalente all’hamiltoniana (1.2.1) del modello di Ising con accoppiamento

Jn,n+µ = |φn · r| · |φn+µ · r| (2.1.4)

Per determinare appieno l’algoritmo bisogna quindi specificare, come abbiamo visto nelcapitolo precedente, i coefficienti Wij del fattore di Boltzmann del modello aumentato

B̂STDn,µ =


)·(W11 W12

W21 W22

)·(ln,µl̄n,µ

). (2.1.5)

11

12 CAPITOLO 2. L’ALGORITMO

Figura 2.1.1: Rappresentazione grafica dei diversi termini del fattore di Boltzmann delmodello aumentato per l’azione Standard.

Figura 2.1.2: Espansione del fattore di Boltzmann del modello aumentato B̂STDσ . Ogni casel-

la colorata corrisponde ad un gruppo di configurazioni. Tutte le configurazioni nello stessogruppo avranno lo stesso peso statistico. Per ogni colonna può essere “attivato” un sologruppo; gli altri avranno peso nullo.

Il diagramma in Figura 2.1.1 mostra il sigificato di ciascuno dei coefficienti Wij . Si vuoleprocede ora in modo diverso rispetto al capitolo precedente per il calcolo dei Wij . Innan-zitutto si nota che la richiesta che per ln,µ = 0 gli spin siano completamente disaccop-piati porta a dire (Equazione 1.2.10) che gli elementi della seconda colonna della matriceW = {Wij} hanno lo stesso peso, che verrà indicato con W2 ed quindi W12 = W22 = W2.Allo stesso modo la condizione che per ln,µ = 1 le orientazioni relative rimangano “con-gelate” (Equazione 1.2.11) equivale ad una esclusione di uno dei due termini della primacolonna di W . Sia W1 il “peso” della prima colonna, a questo punto si hanno solo duepossibilità: W11 = W1 e W21 = 0 oppure W12 = W1 e W11 = 0. Alla luce di questeconsiderazioni la matrice W può essere scritta come il prodotto di una matrice binaria Rche rappresenta le regole secondo cui sono attivati i link al variare dell’orientazione relativadegli spin, ed un “vettore dei pesi”

W =

(W11 W12

W21 W22

)=

(R11 R12

R21 R22

)·(W1

W2

)(2.1.6)

Il diagramma in Figura 2.1.2 mostra una rappresentazione grafica della matrice delle regoleR. Le due diverse opzioni discusse corrispondono ad una diversa scelta del gruppo da“attivare” nella prima colonna.

2.1. AZIONE STANDARD 13

Il problema ora si sposta sulla determinazione dei coefficienti W1 e W2. L’ultimacondizione da imporre consiste nella (1.2.8) che ripetiamo per comodità.∑

ln,µ

B̂STDn,µ = kn,µ(σ)Bn,µ (2.1.7)

Riscriviamo quindi il fattore di Boltzmann come

Bn,µ = p(1)n,µδσn,σn+µ + p(2)n,µδ̄σn,σn+µ (2.1.8)

e notiamo che la condizione (2.1.7) equivale a scrivere(p(1)n,µ

p(2)n,µ

)=

(R11 R12

R21 R22

)·(W1

W2

)= R ·

(W1

W2

)(2.1.9)

e quindi (W1

W2

)= R−1 ·

(p(1)n,µ

p(2)n,µ

)(2.1.10)

Per l’azione standard (2.1.1) abbiamo due possibili regole

R(1) =

(1 10 1

)(2.1.11)

e

R(2) =

(0 11 1

)(2.1.12)

che dati p(1)n,µ = e−βJn,n+µ e p(2)n,µ = e+βJn,n+µriconducono ai casi (1.2.12) e (1.2.13) ottenutiper il modello di Ising.

2.1.1 L’updateA questo punto, data una regolaR valida sotto forma di matrice binaria e dati i coefficientidi accoppiamento p(i)n,µ, l’algoritmo è completamente specificato. L’update della configura-zione del reticolo consiste nella costruzione di un cluster a partire da un sito scelto casual-mente all’interno del reticolo e quindi nel ribaltamento di tutte le variabili σn appartenential cluster. Scendendo nel dettaglio, la procedura per la costruzione del cluster consiste,dato un sito n e per una direzione µ, nell’eseguire un test statistico sulla distribuzione diprobabilità dei link ln,µ condizionata dall’orientazione relativa degli spin σn e σn+µ. Sel’esito è positivo si inserisce il sito n+µ nel cluster, altrimenti non si fa nulla. Si itera cosìil procedimento su tutti i siti aggiunti al cluster finché questo non cessa di espandersi. Aquesto punto tutte le variabili σn del cluster vengono ribaltate (σn → −σn) con probabilità1.

Questa tecnica di update differisce in alcuni dettagli dal caso originale di Swendsen eWang:


Figura 2.1.3: Quattro passi successivi dell’algoritmo per l’hamiltoniana standard di unmodello O(3) (all’equilibrio) su un reticolo 256 × 256 con β = 2.0. La rappresentazioneutilizzata per rappresentare il vettore φn consiste nell’associare ad ogni componente delcampo uno dei tre colori principali della rappresentazione classica RGB.

2.2. AZIONE SYMANZIK 15

Figura 2.2.1: Placchetta di tre siti per l’azione Symanzik.

• in questo caso per ogni step dell’algoritmo si costruisce un solo cluster e non siprocede alla partizione in cluster dell’intero reticolo (come nello Swendsen e Wangoriginale). Questa scelta è stata fatta per ragioni di semplicità.

• Una volta costruito il cluster lo si ribalta con probabilità 1, questo perché non ribaltareun cluster risulterebbe in uno spreco di computazioni dovuto alla costruzione delcluster stesso e non contribuirebbe ad aggiungere statistica.

2.2 Azione SymanzikConsideriamo ora l’hamiltioniana modificata proposta da K. Symanzik

HSYMσ = −

∑n,µ

(4

3φn · φn+µ −

1

12φn · φn+2µ

)(2.2.1)

che con la procedura di embedding proposta da Wolff (2.1.2) diventa

HSYMσ = −

∑n,µ

(4

3Jn,n+µσnσn+µ −

1

12Jn,n+2µσnσn+2µ

). (2.2.2)

Poiché la motivazione di base della scelta dell’azione di Symanzik è quella di utilizzareuna migliore approssimazione della derivata prendendo una combinazione linerare dell’ac-coppiamento con i primi vicini e i secondi vicini, sembra naturale considerare un fattore diBoltzmann di forma

BSYMn,µ = eβ(

43aJn,n+µσnσn+µ+

43(1−a)Jn+µ,n+2µσn+µσn+2µ− 1

12Jn,n+2µσnσn+2µ)

= p(1)n,µδσn,σn+µδσn+µ,σn+2µ + p(2)n,µδσn,σn+µδσn+µ,σn+2µ+

+ p(3)n,µδσn,σn+µδσn+µ,σn+2µ + p(4)n,µδσn,σn+µδσn+µ,σn+2µ (2.2.3)

che per ogni sito n e ogni direzione µ include gli spin σn, σn+µ e σn+2µ. Dato quindi unsito n ed una direzione µ, consideriamo una “placchetta” composta da tre siti in un precisoordine: n, n+ µ, n+ 2µ come in Figura 2.2.1.


Figura 2.2.2: Rappresentazione grafica dei diversi termini del fattore di Boltzmann delmodello aumentato per l’azione Symanzik.

I coefficienti p(i)n,µ saranno dati dalle relazioni

p(1)n,µ = e+43β a Jn,n+µ e+

43β (1−a) Jn+µ,n+2µ e−

112β Jn,n+2µ (2.2.4)

p(2)n,µ = e+43β a Jn,n+µ e−

43β (1−a) Jn+µ,n+2µ e+

112β Jn,n+2µ (2.2.5)

p(3)n,µ = e−43β a Jn,n+µ e−

43β (1−a) Jn+µ,n+2µ e−

112β Jn,n+2µ (2.2.6)

p(4)n,µ = e−43β a Jn,n+µ e+

43β (1−a) Jn+µ,n+2µ e+

112β Jn,n+2µ (2.2.7)

ed inoltre notiamo che ogni termine dell’espansione 2.2.3 corrisponde a due delle 23 = 8possibili configurazioni di spin che differiscono tra loro per un inversione globale.

Come si può vedere dalla Figura2.2.2 per una data configurazione di spin ci sono 5possibili configurazioni dei link, questo ci porta a considerare, per questo modello le matriciR delle regole, introdotte già per l’azione standard, come matrici binarie rettangolari 4× 5.Le regole possibili possono essere determinate imponendo che

1. se due spin sono collegati allora la loro orientazione relativa deve rimanere congelata:sono possibili solo ribaltamenti globali degli spin considerati;

2. se due spin non sono collegati allora devono essere completamente disaccoppiati:tutte le possibili configurazioni dei due spin devono avere lo stesso peso;

3. deve essere vera la relazione (2.1.7) ovvero sommando su tutte le configurazioni pos-sibili delle nuove variabili (link) il fattore di boltzmann B̂SYM

σ deve essere riconduci-bile al modello di partenza (BSYM

σ ).

Imponendo queste condizioni si arriva ad una rappresentazione delle regole possibili comein Figura 2.2.3 dove ogni colonna contiene tutte le possibili configurazioni di spin, tra-lasciando quelle ottenibili da queste per inversione globale, ed ogni riga contiene tutte lepossibili configurazioni dei link. Per poter generare le possibili regole si è proceduto di-videndo le colonne in gruppi di una o più configurazioni. Le condizioni (1) e (2) sonosoddisfatte se

• tutte le configurazioni di un gruppo hanno lo stesso peso,


Figura 2.2.3: Espansione del fattore di Boltzmann del modello aumentato B̂SYMσ . Ogni

casella colorata corrisponde ad un gruppo di configurazioni. Tutte le configurazioni nellostesso gruppo avranno lo stesso peso statistico. Per ogni colonna può essere “attivato” unsolo gruppo; gli altri avranno peso nullo.

• per ogni colonna solo un gruppo può avere peso 1, gli altri dovranno avere peso nullo.

L’ultima condizione è automaticamente soddisfatta imponendo che, seR è una regola possi-bile, allora ogni colonna j (e quindi la j-esima configurazione) dovrà avere un peso statisticoWj definito dalla relazione

p(1)n,µ

p(2)n,µ

p(3)n,µ

p(4)n,µ

= R ·

W1

W2

W3

W4

W5

. (2.2.8)

Abbiamo quindi 32 diverse regole corrispondenti ad una diversa scelta dei gruppi da attivareper ogni colonna. Il risultato è mostrato nella Figura 2.2.4 . Ai fini dell’algoritmo non restaquindi che scegliere una regola ed eseguire l’update.

Sorge ora il problema del calcolo dei coefficienti Wj . Dalla (2.2.8) segue che per ilcalcolo è richiesta l’inversione della matrice R che però è una matrice 4 × 5. I Wj sonoquindi determinati a meno di un parametro λ. Definiamo quindi

R̄ =

(R

0 0 0 0 1

)(2.2.9)

in modo che R̄ sia una matrice quadrata, allora possiamo scrivereW1

W2

W3

W4

W5

= R̄−1 ·

p(1)n,µ

p(2)n,µ

p(3)n,µ

p(4)n,µ

λ

(2.2.10)


A1

B1

C1

D1

E1

F1

G1

H1

A2

B2

C2

D2

E2

F2

G2

H2

A3

B3

C3

D3

E3

F3

G3

H3

A4

B4

C4

D4

E4

F4

G4

H4

Figura 2.2.4: I diagrammi che corrispondono alle regole di update per un algoritmo allaSwendsen e Wang applicato ad un modello con hamiltoniana Symanzik. In rosso sonoraggruppate le regole aventi lo stesso peso e attivabili contemporaneamente. I gruppi consfondo scuro rappresentano le opzioni attive nella regola.


dove λ è un parametro libero non nullo (affinché R̄−1 sia definita). Notiamo che in generaletutti i Wj dipendono da λ ed inoltre, per costruzione W5 = λ.

2.2.1 Studio di un caso particolare: la regola H2

Nell’articolo [4] gli autori espongono un algoritmo particolarmente efficiente, dal punto divista della riduzione del critical slowing-down, per la simulazione di un modello sottopostoall’azione Symanzik. Si vede subito che il caso trattato corrisponde, secondo i metodiappena discussi, ad una particolare scelta di R, a e λ. Si trascuri inizialmente il parametroλ. Vediamo che, l’algoritmo utilizzato nell’articolo [4] corrisponde a scegliere la regolaH2,ovvero

R =

1 1 1 0 01 0 1 1 11 0 0 0 01 1 0 1 0

(2.2.11)

con

a = 1. (2.2.12)

in modo da considerare l’accoppiamento tra primi vicini solo tra σne σn+µ ponendo a zero,nel fattore di Boltzmann l’accoppiamento tra σn+µ e σn+2µ. Per questo, ricordando ilsignificato degli elementi della matrice R (vedi Figura 2.2.3), ci aspettiamo che, al più perun opportuna scelta di λ, si abbia W2 = 0. Vediamo subito che, defniti

Z1 ≡ e−43β Jn,n+µ (2.2.13)

Z2 ≡ e−112β Jn,n+2µ (2.2.14)

si possono scrivere i p(i)n,µ come

p(1)n,µ =Z2

Z1

(2.2.15)

p(2)n,µ =1

Z1Z2

(2.2.16)

p(3)n,µ = Z2Z1 (2.2.17)

p(4)n,µ =Z1

Z2

(2.2.18)


e da queste, attraverso la relazione 2.2.10 si ottiene

W1 = Z1Z2 (2.2.19)

W2 =1

2

[λ− (1− Z2

1) (1− Z22)

Z1Z2

](2.2.20)

W3 =1

2

[(1− Z2

1) (1 + Z22)

Z1Z2

− λ]

(2.2.21)

W4 =1

2

[(1 + Z2

1) (1− Z22)

Z1Z2

− λ]

(2.2.22)

W5 = λ (2.2.23)

e concludiamo che, dovendo essere Wj ≥ 0 per ogni coppia Z1, Z2 possiamo scegliere λnell’intervallo

(1− Z21) (1− Z2

2)

Z1Z2

≤ λ ≤ min

[(1− Z2

1) (1 + Z22)

Z1Z2

,(1 + Z2

1) (1− Z22)

Z1Z2

](2.2.24)

inoltre otteniamo la condizione W2 = 0 se imponiamo che λ sia uguale al suo minimo.Quindi il caso considerato nell’articolo [4] è ottenibile scegliendo:

a = 1 (2.2.25)R = H2 (2.2.26)

λ = λmin =(1− Z2

1) (1− Z22)

Z1Z2

(2.2.27)

2.2.2 Studio delle regole: classi ed applicabilitàSi è visto che le regole che riassumono il metodo di update, sono determinate a meno di unparametro a (il parametro che divide l’interazione tra primi vicini) ed un parametro liberoλ. Le regole inoltre sono in totale 32. Scelta una regola e fissato il parametro a, la richiestache i pesi Wj siano positivi restringe la scelta del parametro λ ad un intervallo. Bisogna oracontrollare se quest’intervallo esiste per ogni valore di a e per ogni regola. L’esistenza omeno di un intervallo valido in cui scegliere λ comporta l’applicabilità o meno della regola.La discussione sull’applicabilità di una regola può essere notevolmente semplificata, anchein modo di velocizzare la computazione, dividendo le regole, in modo opportuno, in classiin modo che data una regola della classe, attraverso opportune trasformazioni è possibilericostruire tutta la classe.

2.2.2.1 Classi di regole e altre semplificazioni

Dalla tabella delle regole in Figura 2.2.4 si nota che il numero di placchette attivate per ogniregola e quindi di 1 nelle rispettive matrici non varia da regola a regola. Questa simmetria


suggerisce che ci si possa ridurre da una matrice ad un altra attraverso semplici permutazionidi righe e colonne in quanto queste lasciano invariati gli elementi delle matrici cambiandolisolo di posto. In particolare, si scambia la riga i-esima e la riga j-esima si riottiene la regoladi partenza a patto di scambiare allo stesso modo l’ordine dei pesi pi e pj . Allo stesso modoscambiando la i-esima colonna e la j-esima colonna si riottiene la stessa regola, a patto discambiare i pesi Wi e Wj . Considerando queste trasformazioni si partiziona l’insieme delle32 regole in 4 classi.

Oltre a queste classi si riesce a riconoscere una struttura più fine considerando la sim-metria sulle regole indotta dalla simmetria per riflessioni del reticolo rispetto ad una delledue direzioni possibli. Si nota che l’inversione di una placchetta corrisponde, in terminidel fattore di Boltzmann, alla trasformazione a→ 1− a ossia ad una inversione rispetto alvalore a = 1/2.

Consideriamo ora l’inversione delle placchette, questa corrisponde ad una particolaretrasformazione delle regole e quindi delle matrici che le rappresentano: si scambiano traloro la seconda e la terza colonna e si scambiano la seconda e la quarta riga. Diciamoquindi che una regola è duale ad un altra se le due regole sono connesse da questa particolaretrasformazione.

Se queste simmetrie sono esatte si dovrà avere che:

1. per ciascuna regola esiste la regola duale a questa (possono anche coincidere)

2. se una regola è applicabile in a allora la sua duale dovrà essere applicabile in 12− a.

In questo modo si semplifica ed in un certo modo si dimezza il numero di controlli da faresulle regole e la loro applicabilità.

Prendiamo ora ciascuna delle quattro classi trovate, determiniamo le trasformazioni trale classi ed il metodo per il calcolo dell’intervallo dei possibili λ.

La prima classe. A questa classe appartengono le regole A1, A2, A4, D1, D3, D4, F1, F2,F3, G2, G3 e G4. Si sceglie la prima regola in ordine alfabetico come rappresentante dellaclasse; le trasformazioni che portano da ciascuna regola al rappresentante della classe sonoriassunte nella Tabella 2.1. Per questa classe si ha quindi

Wj1 =1

2(pi1 + pi2 − pi3 + pi4 − λ)

Wj2 =1

2(−pi1 + pi2 + pi3 − pi4 + λ)

Wj3 =1

2(−pi1 − pi2 + pi3 + p+ λ) (2.2.28)

Wj4 =1

2(pi1 − pi2 + pi3 − pi4 − λ)

Wj5 = λ.


Regole(j1 j2 j3 j4 j5

) (i1 i2 i3 i4

)Regole Duali

A1

(1 2 3 4 5

) (1 2 3 4

)A1

A2

(1 3 4 2 5

) (2 4 3 1

)A4

A4

(1 2 4 3 5

) (4 3 2 1

)A2

D1

(1 2 4 3 5

) (1 3 2 4

)F1

D3

(1 3 4 2 5

) (3 1 2 4

)F3

D4

(1 2 3 4 5

) (4 3 2 1

)F2

F1

(1 3 4 2 5

) (1 3 4 2

)D1

F2

(1 2 3 4 5

) (2 1 4 3

)D4

F3

(1 2 4 3 5

) (3 1 4 2

)D3

G2

(1 2 4 3 5

) (2 4 1 3

)G4

G3

(1 2 3 4 5

) (3 4 1 2

)G3

G4

(1 3 4 2 5

) (4 2 1 3

)G2

Tabella 2.1: Le permutazioni che devono essere utilizzate per riottenere, tramite leEquazioni (2.2.28) e (2.2.29) le regole della prima classe.


) (i1 i2 i3 i4

)Regole Duali

A3

(1 2 3 4 5

) (1 2 3 4

)A3

D2

(1 2 4 3 5

) (1 3 2 4

)F4

F4

(1 3 4 2 5

) (1 3 4 2

)D2

G1

(1 2 4 3 5

) (2 4 1 3

)G1

Tabella 2.2: Le permutazioni che devono essere utilizzate per riottenere, tramite leEquazioni (2.2.30) e (2.2.31) le regole della seconda classe.

La condizione che i pesi siano tutti maggiori di zero porta a restringere le possibili scelte diλ nell’intervallo

max [0, pi1 − pi3 + |pi2 − pi4 |] ≤ λ ≤ pi1 − |pi2 − pi3 + pi4| (2.2.29)

e da questa, utilizzando le permutazioni dei pesi scritte nella tabella si ottengono i Wj el’intervallo di applicabilità per ogni regola della classe.



) (i1 i2 i3 i4

)Regole Duali

B1

(1 2 3 4 5

) (1 2 3 4

)B1

C4

(1 2 3 4 5

) (4 3 2 1

)E2

E2

(1 2 3 4 5

) (2 1 4 3

)C4

H3

(1 2 3 4 5

) (3 4 1 2

)H3

Table 2.3: Le permutazioni che devono essere utilizzate per riottenere, tramite le Equazioni(2.2.33) e (2.2.34) le regole della terza classe.

La seconda classe. La seconda classe è popolata dalle regoleA3, D2, F4 eG1. Prendendocome rappresentante A3 si ottiene

Wj1 =1

2(pi1 + pi2 − pi3 + pi4 + λ)

Wj2 =1

2(−pi1 + pi2 + pi3 − pi4 − λ)

Wj3 =1

2(−pi1 − pi2 + pi3 + p− λ) (2.2.30)

Wj4 =1

2(pi1 − pi2 + pi3 − pi4 − λ)

Wj5 = λ

e

max [0, pi3 − pi1 − pi2 − pi4 ] ≤ λ (2.2.31)pi3 + min [−pi1 + p− pi4 , +pi1 − pi2 − pi4 , −pi1 − pi2 + pi4 ] ≥ λ (2.2.32)

Attraverso le permutazioni in Tabella (2.2) si riottengono tutte le regole della classe.

La terza classe. Alla terza classe appartengono le regole B1, C4, E2 ed H3. Utilizzandola regola B1 otteniamo i seguenti pesi

Wj1 = pi1 − λ

Wj2 =1

2(−pi1 + pi2 + pi3 − pi4 + λ)

Wj3 =1

2(−pi1 − pi2 + pi3 + pi4 + λ) (2.2.33)

Wj4 =1

2(−pi1 + pi2 − pi3 + pi4 + λ)

Wj5 = λ



) (i1 i2 i3 i4

)Regole Duali

B2

(1 2 3 4 5

) (1 2 3 4

)B4

B3

(1 2 4 3 5

) (1 3 2 4

)B3

B4

(1 3 2 4 5

) (1 4 3 2

)B2

C1

(1 3 2 4 5

) (4 1 2 3

)E1

C2

(1 2 4 3 5

) (4 2 3 1

)E4

C3

(1 2 3 4 5

) (4 3 2 1

)E3

E1

(1 2 3 4 5

) (2 1 4 3

)C1

E3

(1 3 2 4 5

) (2 3 4 1

)C3

E4

(1 2 4 3 5

) (2 4 1 3

)C2

H1

(1 2 4 3 5

) (3 1 4 2

)H1

H2

(1 3 2 4 5

) (3 2 1 4

)H4

H4

(1 2 3 4 5

) (3 4 1 2

)H2

Table 2.4: Le permutazioni che devono essere utilizzate per riottenere, tramite le Equazioni(2.2.35) e (2.2.36) le regole della quarta classe.

e da questi l’intervallo dei λ possibili diventa

max (pi1 − pi2 − pi3 + pi4 , pi1 + pi2 − pi3 − pi4 , pi1 − pi2 + pi3 − pi4 , 0) ≤ λ ≤ pi1(2.2.34)

Le trasformazioni che riconducono le Equazioni (2.2.33)(2.2.34) ricavate per la classe B1

alle relazioni equivalenti per tutti gli altri membri della classe sono riassunte nella Tabel-la (2.3).

La quarta classe. Le regoleB2,B3,B4, C1, C2, C3,E1,E3,E4,H1,H2 edH4 formano laquarta ed ultima classe. Scegliando la regola B2 come rappresentante otteniamo i seguentipesi

Wj1 = pi1

Wj2 =1

2(−pi1 + pi2 + pi3 − pi4 − λ)

Wj3 =1

2(−pi1 − pi2 + pi3 + pi4 + λ) (2.2.35)

Wj4 =1

2(−pi1 + pi2 − pi3 + pi4 − λ)

Wj5 = λ

e quindi restringiamo la scelta di λ nell’intervallo

max (pi1 + pi2 − pi3 − pi4 , 0) ≤ λ ≤ pi2 − pi1 − |pi3 − pi4| (2.2.36)

2.3. IMPLEMENTAZIONE DELL’ALGORITMO 25

Le permutazioni che devono essere utilizzate per riottenere, tramite le Equazioni (2.2.35)e (2.2.36) le regole della quarta classe sono riassunte nella Tabella (2.4).

2.3 Implementazione dell’algoritmoL’algoritmo è stato implementato in un software composto essenzialmente da una par-te compilata e da una interpretata: la parte compilata consiste in una libreria dinamica“SWENDSENWANGLIB”, scritta in C++ ed ottimizzata per architetture x86_64 a 64bit; laparte interpretata consiste in alcuni script in PYTHON che utilizzano la libreria SWENDSEN-WANGLIB per eseguire gli update, ottenere e salvare statistica e manipolare dati e risultati.In questo modo si combinano l’efficienza di linguaggi compilati, che permettono numero-se ottimizzazioni in fase di compilazione, e la flessibilità di linguaggi interpretati, come ilPYTHON, utili soprattutto quando si necessita, come in questo caso, di eseguire un numeromolto elevato di processi, con compiti ed impostazioni iniziali anche molto diversi tra loro.

Lo schema che si è seguito è interamente basato sul paradigma della programmazionead oggetti; in particolare si sono definite le seguenti classi:

• SWVECTOR: consente di gestire un vettore di dimensione a scelta e definisce metodiappositi per le operazioni più comuni come la copia, la somma, la moltiplicazione peruno scalare, una rotazione di π e il prodotto scalare tra due vettori;

– SWFIELD: necessaria per la gestione dei siti del reticolo, questa classe è figliadi swVector ed oltre ad ereditare membri e attributi della classe swVector alfine di gestire il vettore dei campi su ogni sito, definisce attributi e metodi perpoter accedere ai primi vicini e per potersi spostare nel reticolo senza bisognodi essere a conoscenza della struttura del reticolo stesso;

• SWLATTICE: definisce il reticolo come un array di siti (swField); la struttura dellattice viene “disegnata” collegando i siti tra loro a piacere. In questa tesi il reticolo èquadrato, ogni sito ha 4 primi vicini, 2 per direzione, e vengono applicate condizioniperiodiche al contorno che fanno sì che il reticolo abbia la topologia di un toro;

• CLUSTERALGORITHM: questa classe definisce metodi e attributi comuni a tutti glialgoritmi che si basano sulla costruzione di cluster (si definiscono qui anche i metodiper acquisire statistiche sui cluster: centroide, dimensione, numero di siti collegati edaltro);

– GENERICSW: con questa classe si definisce il vero e proprio metodo di updateutilizzato, del tutto generale, basato sulla struttura di regole già discussa maindipendentemente dal tipo di azione considerata. Le operazioni e le strutturenecessarie a specializzare l’algoritmo su un modello dato vengono definite nelleclassi figlie di genericSW;


∗ RULESYMSW: questa classe specializza genericSW per un azione Syman-zik e per le regole discusse precedentemente;

• Classi per l’acquisizione di misure e statistica.

Lo schema utilizzato nella costruzione del software è riproposto in Figura 2.3.1 ed in Figura2.3.2.

2.3. IMPLEMENTAZIONE DELL’ALGORITMO 27

Figura 2.3.1: Diagramma delle classi che implementano l’algoritmo di update.


Figura 2.3.2: classiDiagramma delle classi utili all’acquisizione di statistica.

Capitolo 3

Risultati

3.1 Misure da eseguireAl fine di determinare le caratteristiche degli algoritmi e di misurarne l’efficenza, si eseguo-no sui reticoli, sui quali vengono applicati tali algoritmi, tre misure principali [5].

E(1) =1

2

∑n,µ

φn · φn+µ (3.1.1)

M2 =∑n,m

φn · φm (3.1.2)

F =1

2

[∑n,m

e2πı(nx−mx)/Lφn · φm

]+

+1

2

[∑n,m

e2πı(ny−my)/Lφn · φm

](3.1.3)

La prima di queste tre, E(1), rappresenta il primo termine dello sviluppo dell’energia sulreticolo e viene effettuata, non ai fini di una quantificazione del fenomeno del criticalslowing-down ma perché consente, attraverso un confronto diretto con i dati presenti inletteratura, un primo controllo abbastanza pratico e diretto sulla distribuzione di equilibriodell’algoritmo implementato.

Le misure diM2 e F consentono invece di quantificare il fenomeno dello slowing-down.Attraverso queste definiamo

χ = V −1⟨M2

⟩(3.1.4)

χ′ = V −1 〈F 〉 (3.1.5)

29

30 CAPITOLO 3. RISULTATI

dove χ è la suscettività e quindi la lunghezza di correlazione

ξ ≡

(χχ′ − 1

)4 sin2

(πL

)

12

. (3.1.6)

Attraverso quindi la misura di M2, F e i rispettivi tempi di autocorrelazione τint,M2 eτint,F si può determinare l’esponente z che caratterizza il critical slowing-down attraverso larelazione (vedi [3])

τ = ξ(β, L)z φ

[ξ(β, L)

L

](3.1.7)

e quindi determinare l’efficenza della soluzione proposta.

Stima del tempo di autocorrelazione. Uno problema rilevante riscontrato durante l’a-nalisi numerica dei risultati delle simulazioni è la stima del tempo di autocorrelazioneintegrato definito attraverso l’Equazione 1.1.7 che riportiamo qui di seguito:

τint,f =1

2+∞∑t=1

ρff (t). (3.1.8)

Il problema della stima di τint,f per un osservabile f risiede nel fatto che per t → ∞ lafunzione di autocorrelazione tende a zero esponenzialmente e da un punto di vista compu-tazionale e statistico ciò implica che a t molto grandi prevalga il rumore statistico a scapitodel “segnale” (la funzione di autocorrelazione). Questa problematica è stata risolta uti-lizzando il metodo proposto da A.D. Sokal in [14] ovvero introducendo una finestra λ(t)autoregolata e di utilizzare la seguente stima per il tempo di autocorrelazione:

τ̂int,f =1

2+∞∑t=1

λ(t)ρff (t). (3.1.9)

dove, in questo caso, la finestra λ(t) è la funzione a scalino

λ(t) =

{1 per t ≤M

0 altrimenti(3.1.10)

ed M è definito come il più piccolo intero tale che M ≥ c · τ̂int,f (M) dove c è una costantetipicamente compresa tra 4 e 6. In questo modo la varianza di τ̂int,f sarà dell’ordine

Var [τ̂int,f ] ≈2(2M + 1)

nτ̂ 2int,f (3.1.11)

3.2. ANALISI DELL’APPLICABILITÀ DELLE REGOLE. 31

0

20

40

60

80

100

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Acc

etta

nza

Reg

ola

Parametro a

Tutte le regole

A1A2A4D1D3D4F1F2F3G2G3G4A3D2F4G1

B1C4E2H3B2B3B4C1C2C3E1E3E4H1H2H4

Figura 3.2.1: Analisi dell’applicabilità delle regole in Figura 2.2.4 al variare del parametro adi splitting dell’interazione tra primi vicini. La percentuale di successo è riportata sull’asseverticale in funzione del valore di a. Notare la simmetria rispetto ad a = 1

2.


0

20

40

60

80

100

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Acc

etta

nza

Reg

ola

Parametro a

Classe A1

A1A2A4D1D3D4F1F2F3G2G3G4

Figura 3.2.2: Analisi dell’applicabilità delle regole della prima classe. La percentuale disuccesso è riportata sull’asse verticale in funzione del valore di a. Notare la simmetriarispetto ad a = 1

2.


0

20

40

60

80

100

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Acc

etta

nza

Reg

ola

Parametro a

Regola A3

A3D2F4G1

Figura 3.2.3: Analisi dell’applicabilità delle regole della seconda classe. La percentualedi successo è riportata sull’asse verticale in funzione del valore di a. Notare la simmetriarispetto ad a = 1

2.


0

20

40

60

80

100

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Acc

etta

nza

Reg

ola

Parametro a

Regola B1

B1C4E2H3

Figura 3.2.4: Analisi dell’applicabilità delle regole della terza classe. La percentuale disuccesso è riportata sull’asse verticale in funzione del valore di a. Notare la simmetriarispetto ad a = 1

2.


0

20

40

60

80

100

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Acc

etta

nza

Reg

ola

Parametro a

Regola B2

B2B3B4C1C2C3E1E3E4H1H2H4

Figura 3.2.5: Analisi dell’applicabilità delle regole della quarta classe. La percentuale disuccesso è riportata sull’asse verticale in funzione del valore di a. Notare la simmetriarispetto ad a = 1

2.


3.2 Analisi dell’applicabilità delle regole.L’analisi presentata nelle Figure 3.2.1-3.2.5 è stata eseguita testando l’esistenza di un inter-vallo valido per il parametro libero λ al variare del parametro di a di splitting dell’intera-zione tra primi vicini su un numero finito di punti nell’intervallo [−1.5, 2.5]. Le regole sonostate testate su un reticolo di prova 600 × 600 portato all’equilibrio termico con il metododell’Heat Bath. Tra le misure effettuate per diversi valori di a sono stati eseguiti diversicicli dell’algoritmo di update del reticolo in modo da ottenere misure scorrelate tra loro.L’accettanza è data in termini di percentuale di successi di una regola.

Dualità. Analizziamo, per semplicità, le regole della terza classe che sono rappresentatedalla regola B1. Dai risultati presentati in Figura 3.2.4 si nota una palese simmetria rispettoad a = 1

2, in particolare, osservando questi risultati, si osserva che B1 è duale a se stessa,

così come H3 mentre le regole C4 ed E2 sono duali tra loro, ovvero si trasformano l’unanell’altra per riflessioni del parametro a rispetto ad 1

2. La stessa struttura è evidente per tutte

le altre classi (Figure 3.2.2, 3.2.3, 3.2.5) e segue lo schema previsto nella sezione precedenteed espresso nelle Tabelle 2.1, 2.2, 2.3 e 2.4.

Regole da testare. I test sulla accettanza di una regola diminuiscono notevolmente ilnumero di variazioni dell’algoritmo da testare. Come si può osservare dalla Figura 3.2.1per a = 1 solo poche regole sono disponibili per un eventuale algoritmo: si tratta delleregole H2 e D2; allo stesso modo per a = 1

2solo la regola H1è applicabile nel 100% dei

casi. Non viene trattato il caso a = 0 perché le regole applicabili in questo punto sono dualialle regole applicabili in a = 1.

3.3 Energia: confronto con l’Heat Bath.La misura dell’energia sul reticolo offre un controllo veloce sulle distribuzioni di equili-brio dei diversi algoritmi. Dai grafici nelle Figure 3.3.1, 3.3.2 e 3.3.3 si nota l’energiaraggiunta, all’equilibrio, con la nuova soluzione, converge con l’Heat Bath e con l’energiaprevista dalla teoria delle perturbazioni (per indicazioni sul calcolo perturbativo dell’energiavedere [1]).

3.3. ENERGIA: CONFRONTO CON L’HEAT BATH. 37

0.75

0.76

0.77

0.78

0.79

0.8

0 2000 4000 6000 8000 10000

Ene

rgia

t

Confronto Energia

Heat BathCalcolo perturbativo

H1 a=0.5

Figura 3.3.1: Confronto dell’energia su reticoli 100 × 100 a β = 2.0 per l’algoritmo pre-sentato utilizzando la regola H1 con a = 1

2e Heat Bath e l’energia prevista con la teo-

ria delle perturbazioni al terz’ordine. L’energia media con una statistica di 105 campioniall’equilibrio è E(H1)

1 = 0.78749(7) per la regola H1 e EHB1 = 0.78771(3) per l’Heat Bath.


0.75

0.76

0.77

0.78

0.79

0.8

0 2000 4000 6000 8000 10000

Ene

rgia

t

Confronto Energia


H2 a=1.0

Figura 3.3.2: Confronto dell’energia su reticoli 100 × 100 a β = 2.0 per l’algoritmo pre-sentato utilizzando la regola H2 con a = 1 e Heat Bath e L’Energia prevista con la teo-ria delle perturbazioni al terz’ordine. L’energia media con una statistica di 105 campioniall’equilibrio è E(H2)

1 = 0.78781(1) per la regola H2 e EHB1 = 0.78771(3) per l’Heat Bath.

3.3. ENERGIA: CONFRONTO CON L’HEAT BATH. 39

0.75

0.76

0.77

0.78

0.79

0.8

0 2000 4000 6000 8000 10000

Ene

rgia

t

Confronto Energia


D2 a=1.0

Figura 3.3.3: Confronto dell’energia su reticoli 100 × 100 a β = 2.0 per l’algoritmo pre-sentato utilizzando la regola D2 con a = 1 e Heat Bath e l’energia prevista con la teoriadelle perturbazioni al terzo ordine. L’energia media con una statistica di 105 campioniall’equilibrio è E(D2)

1 = 0.78700(1) per la regola D2 e EHB1 = 0.78771(3) per l’Heat Bath.


Figura 3.4.1: Grafico tratto dall’articolo citato in [3]. L’algoritmo utilizzato equivale adapplicare la regola H2 fissando il parametro a = 1.

3.4 Critical Slowing-DownLa valutazione quantitativa del fenomeno del critical slowing-down è stata eseguita per leregole H1 con a = 1

2e D2 con a = 1. Per determinare l’esponente dinamico critico z si è

supposto che sia valida la seguente legge di scala per τint:

τint = ξ(β, L)z φ

[ξ(β, L)

L

](3.4.1)

dove φ è una funzione universale ignota dipendente solo dal rapporto ξ/L. In questo mo-do, misurata la dipendenza di τint dalla lunghezza di correlazione ξ, è possibile determinarel’esponente critico z. Vediamo innanzi tutto dal grafico in Figura 3.4.1 la misura effet-tuata in [3] per un algoritmo equivalente all’applicazione, nel nostro caso, della regola H2

fissando a = 1. In [3], l’esponente critico z è stato determinato utilizzando la (3.4.1) eddato

z ∼ 0.3. (3.4.2)

I grafici in Figura 3.4.2 e in Figura 3.4.3 riportano le misure effettuate per le regole H1

fissando il parametro a = 12

e D2 in a = 1. Le simulazioni sono state svolte utilizzandoreticoli 32× 32, 64× 64 per entrambe le regole e 128× 128 solo per la regola H1.

La Figura 3.4.2 mostra come si dispongono i dati estratti dalla simulazione per z =1.5. Studiando la dipendenza da z si può stimare che l’esponente critico che caratterizza ilcritical slowing-down, varia tra z ∼ 1.3 e z ∼ 1.7 per la regola H1.

3.4. CRITICAL SLOWING-DOWN 41

Figura 3.4.2: Verifica dell’ansatz (3.4.1) per la regola H1. Sono stati utilizzati reticoli conL = 32, L = 64 e L = 128. Il grafico corrisponde al valore di z = 1.5 che rende migliorel’accordo tra i dati per ξ/L > 0.2. I valori con ξ/L < 0.2 sono stati scartati dato che peri reticoli più piccoli il valore di ξ è probabilmente troppo piccolo per trovarsi nel regimeasintotico.

Allo stesso modo dalla Figura 3.4.3 si nota che l’esponente critico è vicino ai valori cheassume per un algoritmo a mosse locali. Entrambe le regole alternative alla H2 peggioranol’efficienza dell’algoritmo.

Per confermare questi risultati si ha bisogno di una statistica più ampia e di retico-li più grandi, per poter essere certi di trovarsi nel regime asintotico caratterizzato dallalegge(3.4.1).

Queste primo studio esplorativo suggerisce che la correlazione tra le probabilità diattivazione dei link (presente nella regola H1) può peggiorare l’efficienza dell’algoritmo.

Lo studio presentato in questa tesi è limitato alle sole tre regole sempre valide in uno deidue valori di a scelti (a = 1 e a = 1

2) e non considera alternative valide con il parametro a

diverso dai valori scelti (vedere la Figura 3.2.1). Un altra variazione possibile sarebbe quelladi considerare una certa combinazione di regole: si applicano prima le regole più efficientie che tendono a generare cluster di dimensioni ridotto: nei casi in cui queste falliscono siutilizza una regola valida.

Tutte queste variazioni sono già previste e implementate nel software sviluppato.


Figura 3.4.3: Verifica dell’ansatz (3.4.1) per la regola D2. Sono stati utilizzati reticoli conL = 32 e L = 64. Il grafico corrisponde al valore di z = 1.85 che rende migliore l’accordotra i dati.

Bibliografia

[1] B. Allés, A. Buonanno, and G. Cella. Perturbation theory predictions and monte carlosimulations for the 2d o(n) non-linear σ-models. Nuclear Physics B, 500(1–3):513 –543, 1997.

[2] K. Binder and D.W. Heermann. Monte Carlo Simulation in Statistical Physics: AnIntroduction. Springer Series in Solid-State Sciences. Springer, 2010.

[3] Alessandra Buonanno and Giancarlo Cella. Swendsen-wang update algorithm for thesymanzik improved σ model. Phys. Rev. D, 51:5865–5869, May 1995.

[4] Alessandra Buonanno, Giancarlo Cella, and Giuseppe Curci. Lattice energy-momentum tensor with symanzik improved actions. Phys. Rev. D, 51:4494–4502,Apr 1995.

[5] Sergio Caracciolo, Robert G. Edwards, Andrea Pelissetto, and Alan D. Sokal. Wol-ff type embedding algorithms for general nonlinear sigma models. Nucl.Phys.,B403:475–541, 1993.

[6] J. Cardy. Scaling and Renormalization in Statistical Physics. Cambridge LectureNotes in Physics. Cambridge University Press, 1996.

[7] Daniel Kandel, Radel Ben-Av, and Eytan Domany. Cluster monte carlo dynamics forthe fully frustrated ising model. Phys. Rev. B, 45:4700–4709, Mar 1992.

[8] J. G. Kemeny and J. L. Snell. Finite Markov Chains. Springer, 1992.

[9] S. Ma. Modern theory of critical phenomena. Number v. 46 in Frontiers in physics.W. A. Benjamin, Advanced Book Program, 1976.

[10] Gene F. Mazenko and Oriol T. Valls. Dynamic critical exponent z in sometwo-dimensional models. Phys. Rev. B, 24:1419–1428, Aug 1981.

[11] Robert B. Pearson, John L. Richardson, and Doug Toussaint. Dynamic correlations inthe three-dimensional ising model. Phys. Rev. B, 31:4472–4475, Apr 1985.

43

44 BIBLIOGRAFIA

[12] M. B. Priestley. Spectral Analysis and Time Series. Academic, 1981.

[13] J. Smit. Introduction to Quantum Fields on a Lattice. Cambridge Lecture Notes inPhysics. Cambridge University Press, 2002.

[14] A.D. Sokal and Troisième cycle de la physique en Suisse romande. Monte CarloMethods in Statistical Mechanics: Foundations and New Algorithms. Troisième Cyclede la Physique. Université de Lausanne, Bâtiment des sciences de physique, 1989.

[15] Robert H. Swendsen and Jian-Sheng Wang. Nonuniversal critical dynamics in montecarlo simulations. Phys. Rev. Lett., 58:86–88, Jan 1987.

Date post:	09-Jul-2020
Category:	Documents
Upload:	others
View:	4 times
Download:	0 times

Un algoritmo di Monte Carlo per la simulazione di un …...un algoritmo di Monte Carlo e riduce...

Documents