基于matlab用蒙特卡罗方法模拟无外加磁场下的二维伊辛模型，绘制能量、磁化强度、热容和磁化率随温度的变化

二维伊辛模型是一个重要的统计物理模型，用于描述强关联体系的基态和激发态。在无外加磁场下，其哈密顿量可以表示为：

$$
H = -J\sum_{<i,j>}s_is_j
$$

其中 $s_i$ 表示在位置 $i$ 的自旋，取值为 $\pm 1$；$<i,j>$ 表示 $i$ 和 $j$ 为相邻格点，$J$ 是自旋之间的交换作用强度。

蒙特卡罗方法是一种常用的计算机算法，可以模拟随机过程，因此非常适合用于二维伊辛模型的模拟。下面我们基于 Matlab 实现无外加磁场下的二维伊辛模型的模拟，并绘制能量、磁化强度、热容和磁化率随温度的变化。

首先，我们需要初始化模型。我们选取一个 $L\times L$ 的正方形格点，每个格点上有一个自旋，初始状态可以随机生成。为了方便计算，我们可以使用周期性边界条件，即相邻的格点相互连接，形成一个环。

% 初始化二维伊辛模型
L = 32; % 格点数
T = linspace(1.5,3,30); % 温度范围
J = 1; % 交换作用强度
N = L*L; % 自旋数
s = sign(rand(N,1)-0.5); % 随机生成初始状态

接下来，我们需要使用 Metropolis 算法进行模拟。具体地，我们每次随机选取一个自旋，计算翻转该自旋后的能量变化 $\Delta E$。如果 $\Delta E<0$，则接受该翻转；否则以一定概率 $\exp(-\Delta E/T)$ 接受该翻转。这里的 $T$ 是温度，控制了模型的热力学行为。

% Metropolis 算法模拟二维伊辛模型
ntherm = 1000; % 热化步数
nmc = 10000; % 蒙特卡罗步数
E = zeros(length(T),nmc); % 能量
M = zeros(length(T),nmc); % 磁化强度
C = zeros(length(T),1); % 热容
X = zeros(length(T),1); % 磁化率
for k=1:length(T)
    for i=1:ntherm
        r = randi([1,N]); % 随机选择一个自旋
        s0 = s;
        s(r) = -s(r); % 翻转该自旋
        deltaE = 2*J*s(r)*(s(mod(r-2,N)+1)+s(mod(r,N)+1)+s(mod(r,L)+r-L)+s(mod(r-1,L)+r-1-L));
        if deltaE > 0 && rand > exp(-deltaE/T(k))
            s(r) = s0(r); % 恢复原状态
        end
    end
    E(k,1) = -J*sum(s.*(s(mod(1:N-1,N)+1)+s(mod(1:N-L-1,N)+L+1)));
    M(k,1) = sum(s)/N;
    for i=2:nmc
        r = randi([1,N]); % 随机选择一个自旋
        s0 = s;
        s(r) = -s(r); % 翻转该自旋
        deltaE = 2*J*s(r)*(s(mod(r-2,N)+1)+s(mod(r,N)+1)+s(mod(r,L)+r-L)+s(mod(r-1,L)+r-1-L));
        if deltaE > 0 && rand > exp(-deltaE/T(k))
            s(r) = s0(r); % 恢复原状态
        end
        E(k,i) = E(k,i-1) + deltaE;
        M(k,i) = M(k,i-1) + 2*s(r)/N;
    end
    E2 = E(k,:).^2;
    M2 = M(k,:).^2;
    C(k) = (mean(E2)-mean(E(k,:))^2)/(T(k)^2);
    X(k) = (mean(M2)-mean(M(k,:))^2)/T(k);
end

最后，我们可以绘制能量、磁化强度、热容和磁化率随温度的变化曲线。这里我们使用 Matlba 的 `plot` 函数进行绘图。

% 绘制能量、磁化强度、热容和磁化率随温度的变化曲线
figure;
subplot(2,2,1);
plot(T,E/N,'.-');
xlabel('Temperature');
ylabel('Energy per spin');
subplot(2,2,2);
plot(T,abs(M),'.-');
xlabel('Temperature');
ylabel('Magnetization per spin');
subplot(2,2,3);
plot(T,C,'.-');
xlabel('Temperature');
ylabel('Specific heat');
subplot(2,2,4);
plot(T,X,'.-');
xlabel('Temperature');
ylabel('Susceptibility');

运行以上代码，即可得到能量、磁化强度、热容和磁化率随温度的变化曲线。