Ising

Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja Marcin (dyskusja | edycje) z dnia 15:30, 27 gru 2010
(różn.) ← poprzednia wersja | przejdź do aktualnej wersji (różn.) | następna wersja → (różn.)

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)

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}\). 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:


  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{\langle m^4\rangle}{3 \langle m^2\rangle^2}\;\]

Proszę narysować wykresy namagnesowania oraz kumulanty U od temperatury

Odnośniki

  • [1] Tobias Preis Homepage (GPU) pdf
  • D. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 2nd edition, Cambridge University Press, Cambridge, 2005.