metropolis.m

用Monte-Carlo方法解决二维Ising模型

% Metropolis算法
function [newsigma,newHsigma]=Metropolis(sigma,Hsigma,beta,preselect)
%调用格式: [newsigma,newHsigma]=Metropolis(sigma,Hsigma,beta,preselect)
%输入:     sigma - 当前状态sigma;
%          Hsigma - 当前状态的能量FuncH(sigma);
%          beta - 参数,即1/(kT) 
%          preselect - 预选方式,1=等可能预选,2=单点变化预选
%输出:     newsigma - 产生的新的状态newsigma
%          newHsigma - 新状态的能量


n=length(sigma);
newsigma=zeros(n,n);

%等可能预选
if preselect==1
    while 1
        for i=1:n
            newsigma(i,:)=unidrnd(2,1,n)*2-3;
        end
        if norm(newsigma-sigma)>0 %sigma~=newsigma
            break;
        end
    end
    newHsigma=FuncH(newsigma);
end


%单点变化预选
if preselect==2
    k=unidrnd(n);
    l=unidrnd(n);
    newsigma=sigma;
    newsigma(k,l)=-sigma(k,l); %随机变化一个sigma中的值   
    %下面这段计算预选态的能量,只计算改变的格点对应的能量的改变.
    %横向能量改变
    if l==n
        deltax=2*(sigma(k,l-1)*sigma(k,l)+sigma(k,1)*sigma(k,l));
    elseif l==1
        deltax=2*(sigma(k,l+1)*sigma(k,l)+sigma(k,n)*sigma(k,1));
    else
        deltax=2*(sigma(k,l-1)*sigma(k,l)+sigma(k,1+1)*sigma(k,l));  
    end
    if k==1
        deltax=2*deltax;
    end
    %纵向能量改变
    if k==n
        deltay=2*(sigma(k-1,l)*sigma(k,l)+sigma(1,l)*sigma(k,l));
    elseif k==1
        deltay=2*(sigma(k+1,l)*sigma(k,l)+sigma(k,l)*sigma(n,l));
    else
        deltay=2*(sigma(k-1,l)*sigma(k,l)+sigma(k+1,l)*sigma(k,l));  
    end
    if l==1
        deltay=2*deltay;
    end
    newHsigma=Hsigma+deltax+deltay;
end

%计算DeltaH(内能的改变)
%这里是可以改进的,因为使用单点法的话只改了一点点

DeltaH=newHsigma-Hsigma;

%如果DeltaH<0,新状态能量低,以概率1接受
%如果DeltaH<0,以一定概率接受
if DeltaH>0
    A=exp(-beta*DeltaH);
    r=rand(1);
    if r>A
        newsigma=sigma;
        newHsigma=Hsigma;
    end
end

end