一维伊辛模型MATLAB

1一维伊辛模型MATLAB源代码clear;clc;%本次模拟中，J=1，k=1kT=0.1:0.1:5;kT_len=length(kT);num=102;Ising=sign(2*rand(1,num)-1);%初始构型N=2*10^5;%------初始化，预分配内存------%E_total=zeros(1,kT_len);E_squ=zeros(1,kT_len);C=zeros(1,kT_len);M_total=zeros(1,kT_len);fori=1:kT_lenE_sum=0;E_squ_sum=0;M_sum=0;forturn=1:NE_old=E_calc(Ising);ri=randi([1,num],1);ifri==1ri=num-1;endifri==numri=2;endIsing(ri)=-Ising(ri);%选择其中一个原子使其自旋态相反E_new=E_calc(Ising);E_diff=E_new-E_old;ifE_diff0&&exp(-E_diff/kT(i))=rand()Ising(ri)=-Ising(ri);%保持原来的构型不变endifturnN/2E=E_calc(Ising);M=sum(Ising)./num;E_sum=E_sum+E;E_squ_sum=E_squ_sum+E^2;M_sum=M_sum+M;end2endE_total(i)=E_sum/(N/2);E_squ(i)=E_squ_sum/(N/2);C(i)=(E_squ(i)-E_total(i).^2)/(kT(i).^2);M_total(i)=M_sum/(N/2);end%-------画图部分--------%subplot(131),plot(kT,E_total,'--*b'),xlabel('kT'),ylabel('能量');title('温度与能量关系图');subplot(132),plot(kT,C,'--*r'),xlabel('kT'),ylabel('热容');title('温度与热容关系图');subplot(133),plot(kT,M_total,'--*k'),xlabel('kT'),ylabel('磁化强度');title('温度与磁化强度关系图');%------函数文件部分，需另外新建函数文件，并保存------%function[E]=E_calc(Ising)%用于计算该构型下的能量len=length(Ising);E=0;fori=2:len-1E=E-Ising(i)*(Ising(i-1)+Ising(i+1))/2;end