3.3.2. MATLAB实现¶
3.3.2.1. Basic metropolis¶
function [C,U] = metropolis(m,n,T)
% 用Metropolis 算法计算 m*n的晶格规模的Ising模型在绝对温度为T时的
% magnetic susceptibility.
k=1;
beta=1/(k*T);
J=1;
%spins=2*round(rand(m,n))-1;
spins=ones(m,n);
N=1e5;
H=0;
Hv=zeros(1,N);
for i=1:m
for j=1:(n-1)
H = H-spins(i,j)*spins(i,j+1);
end
H = H-spins(i,n)*spins(i,1);
end
for j=1:n
for i=1:(m-1)
H = H-spins(i,j)*spins(i+1,j);
end
H = H-spins(m,j)*spins(1,j);
end
x=randi(m,N,1);
y=randi(n,N,1);
r=rand(N,1);
for s=1:N
i=x(s);
j=y(s);
H_new=H;
I=i;
if j==n
J=1;
else
J=j+1;
end
H_new=H_new+2*spins(i,j)*spins(I,J);
I=i;
if j==1
J=n;
else
J=j-1;
end
H_new=H_new+2*spins(i,j)*spins(I,J);
J=j;
if i==m
I=1;
else
I=i+1;
end
H_new=H_new+2*spins(i,j)*spins(I,J);
J=j;
if i==1
I=m;
else
I=i-1;
end
H_new=H_new+2*spins(i,j)*spins(I,J);
dH=H_new-H;
if dH <= 0
H=H_new;
spins(i,j)=-spins(i,j);
else
A=exp(-beta*dH);
if r(s)<=A
H=H_new;
spins(i,j)=-spins(i,j);
end
end
Hv(s)=H;
%draw(spins)
%pause(0.04)
end
H2=0;
H1=0;
for s=(N/2+1):N
H=Hv(s);
H2=H2 + H * H;
H1=H1 + H;
end
H2=2*H2/N;
H1=2*H1/N;
C=beta*beta*(H2-H1*H1)/(m*n);
U=H1/(m*n);
end
function draw(data)
[dimX,dimY]=size(data);
clf
axis([0 dimX 0 dimY])
for x=1:dimX
for y=1:dimY
if data(x,y) == 1
rectangle('Position', [x-1 y-1 1 1], ...
'facecolor','b');
end
end
end
end
3.3.2.2. Swendsen-Wang Algorithm¶
function [C,U] = sw(m,n,T)
% 用Swendsen-Wang 算法计算 m*n的晶格规模的Ising模型在绝对温度为T时的
% magnetic susceptibility.
k=1;%k=8.6173324e-5;
beta=1/(k*T);
J=1;
p=1-exp(-2*beta*J);
spins=ones(m,n);%
%spins=2*round(rand(m,n))-1;
%bond_vertical(i,j)表示spins(i,j)和spins(i+1,j)之间的bond情况，i=m时
%j+1表示1
bond_horizontal = zeros(m,n);
%bond_horizontal(i,j)表示spins(i,j)和spins(i,j+1)之间的bond情况,j=n时
%i+1表示1
bond_vertical = zeros(m,n);
%mark 用于分组的时候标记一个晶格是否被分组
mark=zeros(m,n);
N=1e3;
Hv=zeros(1,N);
for step=1:N
bond_horizontal = zeros(m,n);
bond_vertical = zeros(m,n);
%初始化连接
for i=1:m
for j=1:n-1
if spins(i,j)==spins(i,j+1) && rand()<p
bond_horizontal(i,j)=1;
end
end
if spins(i,n)==spins(i,1) && rand()<p
bond_horizontal(i,n)=1;
end
end
for i=1:m-1
for j=1:n
if spins(i,j)==spins(i+1,j) && rand()<p
bond_vertical(i,j)=1;
end
end
end
for j=1:n
if spins(m,j)==spins(1,j) && rand()<p
bond_vertical(m,j)=1;
end
end
%初始化连接结束，开始分组（cluster）
mark=zeros(m,n);
l=0;
queue=zeros(m*n*2,1);
for i=1:m
for j=1:n
if mark(i,j)==0
flip=2*round(rand())-1;
l = 1;
queue(2*l-1)=i;
queue(2*l)=j;
mark(i,j)=1;
while l ~= 0
x=queue(2*l-1);
y=queue(2*l);
spins(x,y)=spins(x,y)*flip;
l=l-1;
I=x;
if y==n
J=1;
else
J=y+1;
end
if bond_horizontal(x,y)==1 && mark(I,J)==0
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
I=x;
if y==1
J=n;
else
J=y-1;
end
if bond_horizontal(I,J)==1 && mark(I,J)==0
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
J=y;
if x==m
I=1;
else
I=x+1;
end
if bond_vertical(i,j)==1 && mark(I,J)==0
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
J=y;
if x==1
I=m;
else
I=x-1;
end
if bond_vertical(I,J)==1 && mark(I,J)==0
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
end
end
end
end
%draw(spins)
%pause(0.04)
H=energy(spins);
Hv(step)=H;
end
H2=0;%能量二阶矩
H1=0;%能量一阶矩
for s=(N/2+1):N
H=Hv(s);
H2=H2+H*H;
H1=H1+H;
end
H2=2*H2/N;
H1=2*H1/N;
C=beta*beta*(H2-H1*H1)/(m*n);
U=H1/(m*n);
end
function draw(data)
[dimX,dimY]=size(data);
clf
axis([0 dimX 0 dimY])
for x=1:dimX
for y=1:dimY
if data(x,y) == 1
rectangle('Position', [x-1 y-1 1 1], ...
'facecolor','b');
end
end
end
end
function r=energy(data)
r=0;
[dimX,dimY]=size(data);
for x=1:dimX
for y=1:dimY-1
r=r-(data(x,y)*data(x,y+1));
end
r=r-(data(x,dimY)*data(x,1));
end
for y=1:dimY
for x=1:dimX-1
r=r-(data(x,y)*data(x+1,y));
end
r=r-(data(dimX,y)*data(1,y));
end
end
3.3.2.3. The Modification by Wolff¶
function [C,U] = wolff(m,n,T)
% 用Wolff 算法计算 m*n的晶格规模的Ising模型在绝对温度为T时的
% magnetic susceptibility.
k=1;%k=8.6173324e-5;
beta=1/(k*T);
J=1;
p=1-exp(-2*beta*J);
spins=ones(m,n);
%spins=2*round(rand(m,n))-1;
N=1e4;
H=0;
Hv=zeros(1,N);
for i=1:m
for j=1:(n-1)
H = H-spins(i,j)*spins(i,j+1);
end
H = H-spins(i,n)*spins(i,1);
end
for j=1:n
for i=1:(m-1)
H = H-spins(i,j)*spins(i+1,j);
end
H = H-spins(m,j)*spins(1,j);
end
for step=1:N
mark=zeros(m,n);%用于标记一个点是否在集合C中
l=0;
queue=zeros(m*n*2,1);
i=randi(m);
j=randi(n);
l=1;
queue(2*l-1)=i;
queue(2*l)=j;
mark(i,j)=1;
%先生成翻转子集C
while l~=0
i=queue(2*l-1);
j=queue(2*l);
l=l-1;
I=i;
if j==n
J=1;
else
J=j+1;
end
if mark(I,J)==0 && spins(i,j)==spins(I,J) && rand()<p
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
I=i;
if j==1
J=n;
else
J=j-1;
end
if mark(I,J)==0 && spins(i,j)==spins(I,J) && rand()<p
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
J=j;
if i==m
I=1;
else
I=i+1;
end
if mark(I,J)==0 && spins(i,j)==spins(I,J) && rand()<p
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
J=j;
if i==1
I=m;
else
I=i-1;
end
if mark(I,J)==0 && spins(i,j)==spins(I,J) && rand()<p
l=l+1;
queue(2*l-1)=I;
queue(2*l)=J;
mark(I,J)=1;
end
spins(i,j)=spins(i,j)*(-1);
end
%下面修正能量。
for i=1:m
for j=1:n
if mark(i,j)==1
I=i;
if j==n
J=1;
else
J=j+1;
end
if mark(I,J)==0
H=H-2*spins(i,j)*spins(I,J);
end
I=i;
if j==1
J=n;
else
J=j-1;
end
if mark(I,J)==0
H=H-2*spins(i,j)*spins(I,J);
end
J=j;
if i==m
I=1;
else
I=i+1;
end
if mark(I,J)==0
H=H-2*spins(i,j)*spins(I,J);
end
J=j;
if i==1
I=m;
else
I=i-1;
end
if mark(I,J)==0
H=H-2*spins(i,j)*spins(I,J);
end
end
end
end
Hv(step)=H;
%draw(spins)
%bpause(0.1)
end
H2=0;
H1=0;
for i=(N/2+1):N
H=Hv(i);
H2=H2+H*H;
H1=H1+H;
end
H2=2*H2/N;
H1=2*H1/N;
C=beta*beta*(H2-H1*H1);
C=C/(m*n);
U=H1/(m*n);
end
function draw(data)
[dimX,dimY]=size(data);
clf
axis([0 dimX 0 dimY])
for x=1:dimX
for y=1:dimY
if data(x,y) == 1
rectangle('Position', [x-1 y-1 1 1], ...
'facecolor','b');
end
end
end
end