Z Skrypty dla studentów Ekonofizyki UPGOW
m (→Fast Ising (matlab)) |
m (→Fast Ising (matlab)) |
||
| Linia 91: | Linia 91: | ||
Bierzemy spin w wybranym punkcie <math>s=1</math>, w takim wypadku mamy | Bierzemy spin w wybranym punkcie <math>s=1</math>, w takim wypadku mamy | ||
| - | <math>E_k = - \sum_{i=1,2,3,4} s_k s_i =- \sum_{i=1,2,3,4} s_i | + | <math>E_k = - \sum_{i=1,2,3,4} s_k s_i =- \sum_{i=1,2,3,4} s_i \,</math> |
| - | </math> | + | |
| + | odpowiada to liniom kodu: | ||
| + | <source lang="matlab"> | ||
| + | el = t(e-1) + t(e+1) + t(e-p-2) + t(e+p+2); | ||
| + | t(e) = 1; | ||
| + | eel = ex(5+el); | ||
| + | </source> | ||
| + | które obliczają to równolegle dla całej podsieci parzystej (e). | ||
| + | <source lang="matlab"> | ||
| + | t(e(find(rand(1,nr)<(eel./(eel+1./eel))))) = -1; | ||
| + | </source> | ||
| + | |||
| + | <source lang="matlab"> | ||
| + | t(1,:)=t(p+1,:); t(p+2,:)=t(2,:); t(:,1)=t(:,n+1); t(:,n+2)=t(:,2); | ||
| + | </source> | ||
Całe źródło: | Całe źródło: | ||
Wersja z 10:04, 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)
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}\)
Bierzemy spin w wybranym punkcie \(s=1\), w takim wypadku mamy
\(E_k = - \sum_{i=1,2,3,4} s_k s_i =- \sum_{i=1,2,3,4} s_i \,\)
odpowiada to liniom kodu:
el = t(e-1) + t(e+1) + t(e-p-2) + t(e+p+2); t(e) = 1; eel = ex(5+el);
które obliczają to równolegle dla całej podsieci parzystej (e).
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);
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;