Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 16:17, 16 gru 2010

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)

Paul Fieguth

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}$! 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 zsumować te procesy i przyjać, że dany spin ma wartość 1 wylosować z prawdopodobieństem $P_{+1\to -1}$ jego zmianę na spin -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:

1. 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.
2. 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.
3. Dla poszczególnych temperatur wyznaczyć kumulantę Bindera.

$$U =1 - \frac{<m^4>}{3 <m^2>^2}\;$$

Proszę narysować wykresy namagnesowania oraz kumulanty U od temperatury

Odnośniki

• [1] Tobias Preis Homepage (GPU) a w szczególności prace
• D. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 2nd edition, Cambridge University Press, Cambridge, 2005.