伊辛(ising)模型MATLAB源程序

de=2.*spin(a(1),a(2)).*ds; if rand<exp(-de/T) spin(a(1),a(2))=-spin(a(1),a(2)); Ett=Ett+de; end Mt(i)=sum(sum(spin))/(N^2); Et0(i)=Ett; if mod(i,step)==0 set(gspin,'Cdata',double(spin)); set(gMt,'Ydata',Mt); set(gEt,'Ydata',Et0); drawnow end end Cv(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); drawnow set(gCv,'Ydata',Cv); drawnow set(gxx,'Ydata',xx); drawnow L=L+1; end
clear all; close all; clc N=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; for T=Te Mt=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([0 imax -1.1 1.1]);
Ett=Et0(1); title('固定温度下的能量'); %xlabel('时间'); %ylabel('Et'); %energy end subplot(3,2,4); gMT=plot(Te,MT,'o-'); axis([0 max(Te) 0 1.1]); title('磁矩随温度T的变化'); %xlabel('T'); %ylabel('MT'); subplot(3,2,5) gCv=plot(Te,Cv,'o-'); %axis([0 max(Te) 0 700]); title('热容随温度T的变化'); %xlabel('T'); %ylabel('Cv'); subplot(3,2,6) gxx=plot(Te,xx,'o-'); %axis([0 max(Te) 0 0.1]); title('磁导率随温度T的变化'); %xlabel=('T'); %ylabel=('X');
Baidu Nhomakorabea
for i=2:imax a=ceil(rand(1,2).*N); if a(1)==1 xl=N; else xl=a(1)-1; end if a(1)==N xr=1; else xr=a(1)+1; end if a(2)==1 yd=N; else yd=a(2)-1; end if a(2)==N yu=1; else yu=a(2)+1; end ds=spin(a(1),yu)+spin(a(1),yd)+spin(xl,a(2))+spin(xr,a(2));
title('固定温度下的磁矩'); %xlabel('时间'); %ylabel('Mt'); %initial energy idx=repmat({':'},ndims(spin),1); idx{1}=[N 1:N-1]; Etotal=spin(idx{:}); idx=repmat({':'},ndims(spin),1); idx{2}=[N 1:N-1]; Etotal=Etotal+spin(idx{:}); idx=repmat({':'},ndims(spin),1); idx{1}=[2:N 1]; Etotal=Etotal+spin(idx{:}); idx=repmat({':'},ndims(spin),1); idx{2}=[2:N 1]; 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);