Z Skrypty dla studentów Ekonofizyki UPGOW
Spis treści |
Model Isinga w ekonofizyce
Mamy Hamiltonian w ogólnej postaci: \[H = -\sum_{i,j} J_{ij} s_i s_j - \mu B \sum_i s_i \, \]
Model Isinga jest szczegolnym przypadkiem, w którym spiny \(s_i\) znajdują się na sieci a oddziaływanie jest ograniczone do najbliższych sąsiadów. O ile w ogólnym przypadku problem ma złożoność obliczeniową \(N^2\) to w modelu Isinga ta złożoność jest rzędu \(N\).
Konkretna realizacja jest dana przez zadanie sieci oraz siły oddziaływania \(J\). Wyróżniamy:
- J > 0 oddziaływanie ferromagnetyczne;stan uporządkowany ma najniższą energię
- J < 0 oddziaływanie antyferromagnetyczne;stan o spinach naprzemiennych ma najniższą energię
Fast Ising (matlab)
Opis działania:
Program implementuje metodę symulacji modelu Isinga stosując tzw. checkerboard decomposition polegającą na podziale sieci na dwie niezależne i aktualizacji wszystkich węzłów każdej podsieci w jednocześnie.
- załóżmy, że przeprowadzamy symulacje dla dwuwymiarowego modelu przy n=6
- definiowane są dwie tablice 8x8 (dodajemy po jednym rzędzie spinów z każdej strony):
t = 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 e = find(t);
t = 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 o = find(t);
dla których miejsca (indeksy) niezerowych elementów sa zapamiętanie w wektorach e i o. Proszę zauważyć, że e i o są jednowymiarowymi tablicami wskaźników. W Octave/Matlab macierz można indeksować zarówno jednym jak i dwoma wskaźnikami np.
octave:28> M=[1,2;3,4] M = 1 2 3 4 octave:29> M(1) ans = 1 octave:30> M(1,2) ans = 2 octave:31> M(3) ans = 2 octave:32> M(1:4) ans = 1 3 2 4
Tak więc w wektorach e i o są miejsca podsieci parzystej i nieparzystej w reprezentacji jednowskaźnikowej.
Funkcja exp jest tablicowana:
ex = exp(-b*[-4:4]);
Zadajemy losowy warunek początkowy:
t = sign(randn(p+2,n+2));
Wartości spinów na kazdej z podsieci uaktualniamy stosująć dynamikę Glaubera, w której prawdopodobieństwo przejścia ze stanu \(i\to j\) wynosi:
\[P_{i\to j} = \displaystyle\frac{e^{-E_j \beta}}{e^{-E_j\beta}+e^{-E_i\beta}}\]
Rozważamy stan energetyczny spinu k i jego najbliższych sąsiadów:
\[E_k = - J \sum_{i=1,2,3,4} s_k s_i =- J s_k \sum_{i=1,2,3,4} s_i = -J s_k E_s\,\]
czyli dla spinu z jego sąsiadami mamy (kładąc: \(b=J \beta\,\)):
\[P_{+1\to -1} = \displaystyle\frac{e^{-E_s b}}{e^{-E_s b}+e^{+E_s b}}=\frac{e^{-E_s b}}{e^{-E_s b}+1/e^{-E_s b}}= \frac{eel}{eel+1/eel}\],
gdzie \(eel = e^{-E_s b}\,\).
odpowiada to linii kodu:
el = t(e-1) + t(e+1) + t(e-p-2) + t(e+p+2);
podobnie mamy \[P_{-1\to +1} = \displaystyle\frac{e^{E_s b}}{e^{E_s b}+e^{-E_s b}}= \frac{1/eel}{eel+1/eel}\].
Zauważmy, że zachodzi \(P_{+1\to -1}=1-P_{-1\to +1}\). W konsekwencji prawdopodobienstwo całkowite, że spin jest w stanie +1 wynosi: \[P_{-1\to +1} P_{-1} +(1-P_{+1\to -1}) P_{+1}= P_{-1\to +1} P_{-1}+P_{-1\to +1} P_{+1}= P_{-1\to +1} \] i podobnie dla stanu -1: \[P_{+1\to -1} P_{+1} +(1-P_{-1\to +1}) P_{-1}= P_{+1\to -1} P_{-1}+P_{+1\to -1} P_{+1}= P_{+1\to -1}=1-P_{-1\to +1} \]
Oznacza to, że zamiast sprawdzać jaki jest spin w węźle i odpowiednio stosować \(P_{-1\to +1}\) lub \(P_{+1\to -1}\), możemy wylosować stan, że dany spin ma wartość +1 z prawdopodobieństem \(P_{-1\to +_1}\).
Przekłada się to na następującą implementację:
t(e) = 1; eel = ex(5+el); t(e(find(rand(1,nr)<(eel./(eel+1./eel))))) = -1;
która oblicza to równolegle dla całej podsieci parzystej (e). Z symulacji zakładamy okresowe warunki brzegowe:
t(1,:)=t(p+1,:); t(p+2,:)=t(2,:); t(:,1)=t(:,n+1); t(:,n+2)=t(:,2);
Komentarza wymaga jeszcze rola zmiennej q. Jest to liczba niezależnych wykonać symulacji (każda z oddzielnym warunkiem początkowym), tak by macierz zwracana przez fukcję ising2 zawiera wszystkie te wyniki i ma wymiar (p*q,n).
s = zeros(p*q,n);
Całe źródło:
function s = ising2(n,m,b,q) n=round(n); m=round(m); if (nargin < 4), q=1; else, q=round(q); end; p=n; if (length(n)>1), p=n(1); n=n(2); end; s = zeros(p*q,n); t = zeros(p+2,n+2); for i=2:(p+1), for j=(2+rem(i,2)):2:(n+1), t(i,j)=1; end; end; e = find(t); t = zeros(p+2,n+2); for i=2:(p+1), for j=(3-rem(i,2)):2:(n+1), t(i,j)=1; end; end; o = find(t); nr = length(e); ex = exp(-b*[-4:4]); for qi=1:q, t = sign(randn(p+2,n+2)); for j=1:m, el = t(e-1) + t(e+1) + t(e-p-2) + t(e+p+2); t(e) = 1; eel = ex(5+el); t(e(find(rand(1,nr)<(eel./(eel+1./eel))))) = -1; t(1,:)=t(p+1,:); t(p+2,:)=t(2,:); t(:,1)=t(:,n+1); t(:,n+2)=t(:,2); el = t(o-1) + t(o+1) + t(o-p-2) + t(o+p+2); t(o) = 1; eel = ex(5+el); t(o(find(rand(1,nr)<(eel./(eel+1./eel))))) = -1; t(1,:)=t(p+1,:); t(p+2,:)=t(2,:); t(:,1)=t(:,n+1); t(:,n+2)=t(:,2); end; s((1:p)+(qi-1)*p,:) = t(2:(p+1),2:(n+1)); end;
Zadania
1) Wyznaczanie temperatury krytycznej modelu Isinga (2D) w zerowym polu.
Proszę napisać program, który wyznaczy temperaturę przejścia fazowego. Dla dwuwymiarowego modelu Isinga na sieci kwadratowej rozwiązanie dokładne jest znane i wynosi \[T_c = \frac{2}{\log(1+\sqrt{2})}\] Aby wyliczyć temperaturę krytyczną z sumylacji Monte-Carlo należy:
- Dla wybranych temperatur należy, startując z losowego stanu początkowego ztermalizować układ tj. wysymulwać na tylę dużą liczbę kroków (np. 1e4), by układ był w równowadze.
- Wykonać M kroków, dla których obliczamy namagnetyzowanie m oraz m^2 i m^4. Należy uśrednić wielkości otrzymane w poszczególnuch krokach.
- Dla poszczególnych temperatur wyznaczyć kumulantę Bindera.
\[U =1 - \frac{\langle m^4\rangle}{3 \langle m^2\rangle^2}\;\]
Proszę narysować wykresy namagnesowania oraz kumulanty U od temperatury