Z Skrypty dla studentów Ekonofizyki UPGOW
Model Isinga w ekonofizyce
Mamy Hamiltonian:
Model Isinga dany jest przez ogólną
- ferromagnetic coupling
- anti-ferromagnetic coupling
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{exp(-E_j*\beta)}{exp(-E_j*\beta+exp(-E_i*\beta}
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;