基于matlab用蒙特卡罗方法模拟无外加磁场下的二维伊辛模型

伊辛模型是一个经典的统计物理模型，描述了一系列自旋在二维晶格上的相互作用。在无外加磁场的情况下，伊辛模型的哈密顿量可以表示为：

H = -J ∑<i,j> σiσj

其中，J 是交换耦合常数，<i,j> 表示 i,j 为相邻格点，σi 取值为 ±1，表示在 i 格点处的自旋方向。

蒙特卡罗方法是一种基于随机数的数值计算方法，可以用来模拟复杂的物理系统。在模拟伊辛模型时，蒙特卡罗方法通常采用 Metropolis 算法，其基本思想是通过随机扰动自旋状态，计算新状态与原状态的能量差 ΔE，然后根据 Metropolis 准则接受或拒绝新状态。

下面是一个基于 Matlab 的二维伊辛模型模拟程序：

clear all;

% 设置参数
L = 50;                   % 模拟晶格大小
J = 1;                    % 交换耦合常数
T = 2.0;                  % 模拟温度
N = L*L;                  % 自旋总数
steps = 10000;            % 模拟步数
w = zeros(1,steps);       % 记录自旋翻转概率

% 初始化自旋状态
sigma = ones(L,L);
for i = 1:L
    for j = 1:L
        if rand() < 0.5
            sigma(i,j) = -1;
        end
    end
end

% 开始模拟
E = -J*sum(sum(sigma.*circshift(sigma,[1,0])))-J*sum(sum(sigma.*circshift(sigma,[0,1])));
for k = 1:steps
    % 随机选取一个自旋
    i = randi(L);
    j = randi(L);
    % 计算能量差
    DeltaE = 2*J*sigma(i,j)*(sigma(mod(i,L)+1,j)+sigma(i,mod(j,L)+1)+sigma(mod(i-2,L)+1,j)+sigma(i,mod(j-2,L)+1));
    % 计算自旋翻转概率
    w(k) = exp(-DeltaE/T);
    % Metropolis 准则
    if rand() < w(k)
        sigma(i,j) = -sigma(i,j);
        E = E + DeltaE;
    end
end

% 绘制自旋状态图和自旋翻转概率变化图
figure;
subplot(1,2,1);
imagesc(sigma);
title(['T=',num2str(T),',E=',num2str(E)]);
axis square;
colormap(gray);
subplot(1,2,2);
plot(1:steps,w);
title('Monte Carlo Steps');
xlabel('Step');
ylabel('Acceptance Probability');

在程序中，首先设置模拟的晶格大小 L、交换耦合常数 J、模拟温度 T 等参数。然后，通过随机函数初始化自旋状态 sigma。接着，在模拟过程中，按照 Metropolis 算法的流程，随机选取一个自旋，计算能量差 DeltaE 和自旋翻转概率 w，然后根据 Metropolis 准则接受或拒绝新状态。最后，绘制自旋状态图和自旋翻转概率变化图。

需要注意的是，在 Metropolis 准则的计算中，需要使用 exp 函数来计算自旋翻转概率。在温度较低的情况下，自旋翻转概率很小，容易出现数值下溢的情况，因此需要采用一些技巧来避免这种情况的发生。例如，可以将自旋翻转概率 w 限制在 1 以下，然后使用 rand 函数来模拟随机数的生成。