xiangbian.m

Ising模型铁磁—顺磁相变现象的Monte Carlo模拟的matlab程序

clear
clc
% Ising Model 初始化各磁畴方向
Ising=zeros(22,22);
for i=1:22
    for j=1:22
        if rem(j,2)==0
           Ising(i,j)=-1;
           plot(i,j,'rv')
           hold on
        else 
           Ising(i,j)=1;
           plot(i,j,'k^')
           hold on
        end
     end
end
axis off
% Monte Carlo 重要抽样方法使铁磁处于平衡状态
Tup=0.1:0.1:5.0;
Tdown=5.0:-0.1:0.1;
%T增大
for i=1:50
    E1=0;
    E2=0;
    E3=0;
    M1=0;
    M2=0;
    M3=0;
    for k=1:1000000
        sai=round(20*rand+1.5);
        saj=round(20*rand+1.5);
        Ising(1,:)=Ising(21,:);
        Ising(22,:)=Ising(2,:);
        Ising(:,1)=Ising(:,21);
        Ising(:,22)=Ising(:,2);
        U1=-Ising(sai,saj)*(Ising(sai-1,saj)+Ising(sai+1,saj)
+Ising(sai,saj-1)+Ising(sai,saj+1))/Tup(i);
        U2=-U1;
        if exp(U1-U2)>rand
           Ising(sai,saj)=-Ising(sai,saj);
        end
        if k>500000       %求出T温度下的Ev,Cv,M,X
           Ising(1,:)=Ising(21,:);
           Ising(22,:)=Ising(2,:);
           Ising(:,1)=Ising(:,21);
           Ising(:,22)=Ising(:,2);
           E=0; 
           Mo=0;
           for i1=2:21
               for j1=2:21
                   E=E-Ising(i1,j1)*(Ising(i1-1,j1)+Ising(i1+1,j1)
+Ising(i1,j1-1)+Ising(i1,j1+1));
                   Mo=Mo+Ising(i1,j1);
               end
           end
           E=E/2;
           Mo=Mo/400;
           E1=E1+E;
           E2=E2+1; 
           E3=E3+E^2;
           M1=M1+Mo;
           M2=M2+1;
           M3=M3+Mo^2;
        end
    end
    Ev(i)=E1/E2;
    E2v(i)=E3/E2;
    Cv(i)=(E2v(i)-Ev(i)^2)/(Tup(i)^2);
    M(i)=M1/M2;
    M2v(i)=M3/M2;
    X(i)=(M2v(i)-M(i)^2)/Tup(i);
end
% T减小
for i=1:50
    E11=0;
    E21=0;
    E31=0;
    M11=0;
    M21=0;
    M31=0;
    for k=1:1000000
        sai=round(20*rand+1.5);
        saj=round(20*rand+1.5);
        Ising(1,:)=Ising(21,:);
        Ising(22,:)=Ising(2,:);
        Ising(:,1)=Ising(:,21);
        Ising(:,22)=Ising(:,2);
        U11=-Ising(sai,saj)*(Ising(sai-1,saj)+Ising(sai+1,saj)
+Ising(sai,saj-1)+Ising(sai,saj+1))/Tdown(i);
        U21=-U11;
        if exp(U11-U21)>rand
           Ising(sai,saj)=-Ising(sai,saj);
        end
        if k>500000       %求出T温度下的Ev,Cv,M,X
           Ising(1,:)=Ising(21,:);
           Ising(22,:)=Ising(2,:);
           Ising(:,1)=Ising(:,21);
           Ising(:,22)=Ising(:,2);
           E0=0;
           Moo=0;
           for i1=2:21
               for j1=2:21
                   E0=E0-Ising(i1,j1)*(Ising(i1-1,j1)+Ising(i1+1,j1)
+Ising(i1,j1-1)+Ising(i1,j1+1));
                   Moo=Moo+Ising(i1,j1);
               end
           end
           E0=E0/2;
           Moo=Moo/400;
           E11=E11+E0;
           E21=E21+1; 
           E31=E31+E0^2;
           M11=M11+Moo;
           M21=M21+1; 
           M31=M31+Moo^2;
        end
    end
    Ev0(i)=E11/E21;
    E2v0(i)=E31/E21;
    Cv0(i)=(E2v0(i)-Ev0(i)^2)/(Tdown(i)^2);
    M0(i)=M11/M21;
    M2v0(i)=M31/M21;
    X0(i)=(M2v0(i)-M0(i)^2)/Tdown(i);
end
% 分别作Ev-T,Cv-T,M-T,X-T的关系图
figure
Tup=0.1:0.1:5.0;
Tdown=5.0:-0.1:0.1;
plot(Tup,Ev,'--*k',Tdown,Ev0,'--xr')
xlabel('温度 T');
ylabel('平均能量 Ev');
legend('Tup','Tdown');
title('温度与平均能量的关系图');
figure
plot(Tup,Cv,'--*k',Tdown,Cv0,'--xr')
xlabel('温度 T');
ylabel('热容 Cv');
legend('Tup','Tdown');
title('温度与热容的关系图');
figure
plot(Tup,M,'--*k',Tdown,M0,'--xr')
xlabel('温度 T');
ylabel('磁化强度 M');
legend('Tup','Tdown');
title('温度与磁化强度的关系图');
figure
plot(Tup,X,'--*k',Tdown,X0,'--xr')
xlabel('温度 T');
ylabel('磁化率 X');
legend('Tup','Tdown');
title('温度与磁化率的关系图');