Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 10:33, 14 gru 2010

Model Isinga w ekonofizyce

Mamy Hamiltonian: \(H = -\sum_{i,j} J_{ij} s_i s_j - \mu B \sum_i s_i \, \)

Model Isinga dany jest przez ogólną

• ferromagnetic coupling
• anti-ferromagnetic coupling

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{exp(-E_j \beta)}{exp(-E_j\beta+exp(-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{exp(-E_s b)}{exp(-E_s b+exp(+E_s b}=\frac{exp(-E_s b)}{exp(-E_s b+1/exp(-E_s b}= \frac{eel)}{eel+1/eel}\),

gdzie \(eel = exp(-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{exp(E_s b)}{exp(E_s b+exp(-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);

.

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;

Problem czasu?