伊辛(ising)模型MATLAB源程序

clearall;closeall;clcN=16;Te=[0:0.2:5];imax=500000;MT=zeros(1,length(Te));Cv=zeros(1,length(Te));xx=zeros(1,length(Te));L=1;k=1;forT=TeMt=zeros(1,imax);step=100000;spin=round(rand(N)).*2-1;subplot(3,2,1)gspin=pcolor(spin);colormap(gray(2));title('自旋分布情况');Mt=zeros(1,imax);subplot(3,2,2)gMt=plot([1:imax],Mt);axis([0imax-1.11.1]);title('固定温度下的磁矩');%xlabel('时间');%ylabel('Mt');%initialenergyidx=repmat({':'},ndims(spin),1);idx{1}=[N1:N-1];Etotal=spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{2}=[N1:N-1];Etotal=Etotal+spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{1}=[2:N1];Etotal=Etotal+spin(idx{:});idx=repmat({':'},ndims(spin),1);idx{2}=[2:N1];Etotal=Etotal+spin(idx{:});Etotal=spin.*Etotal;Et0=zeros(1,imax);Et0(1)=-sum(sum(Etotal))/2;subplot(3,2,3);gEt=plot([1:imax],Et0);Ett=Et0(1);title('固定温度下的能量');%xlabel('时间');%ylabel('Et');%energyendsubplot(3,2,4);gMT=plot(Te,MT,'o-');axis([0max(Te)01.1]);title('磁矩随温度T的变化');%xlabel('T');%ylabel('MT');subplot(3,2,5)gCv=plot(Te,Cv,'o-');%axis([0max(Te)0700]);title('热容随温度T的变化');%xlabel('T');%ylabel('Cv');subplot(3,2,6)gxx=plot(Te,xx,'o-');%axis([0max(Te)00.1]);title('磁导率随温度T的变化');%xlabel=('T');%ylabel=('X');fori=2:imaxa=ceil(rand(1,2).*N);ifa(1)==1xl=N;elsexl=a(1)-1;endifa(1)==Nxr=1;elsexr=a(1)+1;endifa(2)==1yd=N;elseyd=a(2)-1;endifa(2)==Nyu=1;elseyu=a(2)+1;endds=spin(a(1),yu)+spin(a(1),yd)+spin(xl,a(2))+spin(xr,a(2));de=2.*spin(a(1),a(2)).*ds;ifrandexp(-de/T)spin(a(1),a(2))=-spin(a(1),a(2));Ett=Ett+de;endMt(i)=sum(sum(spin))/(N^2);Et0(i)=Ett;ifmod(i,step)==0set(gspin,'Cdata',double(spin));set(gMt,'Ydata',Mt);set(gEt,'Ydata',Et0);drawnowendendCv(L)=1/(T^2).*(mean(Et0.^2)-mean(Et0)^2);xx(L)=1/T.*(mean(Mt.^2)-mean(Mt)^2);MT(L)=mean(abs(Mt([0.9*imax:imax])));set(gMT,'Ydata',MT);drawnowset(gCv,'Ydata',Cv);drawnowset(gxx,'Ydata',xx);drawnowL=L+1;end