MATLAB二维Ising模型Monte-Carlo模拟：磁化曲线不收敛的避坑指南

简介：这份资源面向凝聚态物理、材料科学方向的学习者与科研人员，提供一套用MATLAB编写的二维Ising模型Monte-Carlo模拟与磁化特性分析程序，帮助在无需复杂物理实验的条件下研究铁磁材料的相变行为。压缩包共2个文件，约4KB，包含1个m脚本与1个md说明文档：m文件承载核心模拟逻辑，实现二维自旋网格初始化、Metropolis抽样迭代以及磁化强度与能量的统计计算；md文档则交代项目背景与使用说明，便于快速上手。程序可输出平均磁化强度和能量随温度变化的曲线，并支持可视化不同温度下自旋网格的排列状态，直观呈现临界温度附近的显著变化。目前已有105人学习下载，适合希望理解Ising模型、Metropolis算法与统计物理数值模拟的读者参考，也可作为相关课程设计或科研入门的工具脚本。

1. 从一次“磁化曲线不收敛”说起：二维 Ising 模型到底在算什么

很多人第一次在 MATLAB 里跑二维 Ising 模型，都会遇到一个诡异现象：温度设到 2.0 附近，磁化强度曲线像喝醉了酒，一会儿冲到 0.9，一会儿掉到 0.1，换一组随机数种子结果又完全不一样。这不是代码写错了，而是 Monte-Carlo 模拟本身就在跟临界慢化、有限尺寸效应和采样自相关较劲。二维 Ising 模型是统计物理里最经典的格点自旋系统，每个格点取 ±1，只和最近邻耦合，却能自发涌现出二级相变，Onsager 在 1944 年就给出了严格解，临界温度约 2.269。用 MATLAB 实现基于 Monte-Carlo 模拟的二维 Ising 模型磁化分析系统，核心目标不是复现解析解，而是搭一套能跑出磁化强度、能量、比热、磁化率随温度变化曲线的可复现流程，并理解 Metropolis 采样为什么在临界点附近会“翻车”。这套东西适合谁？做计算物理课程设计的学生、想用 MATLAB 验证相变理论的工程师、以及需要把随机模拟结果和解析解对齐的从业者。下面我按自己踩过的坑，把选型、实现、参数和排查一条条讲清楚。

2. 用 MATLAB 搭 Metropolis 采样：从单自旋翻转写到整条磁化曲线

2.1 为什么选 Metropolis 而不是 Wolff 或 Swendsen-Wang

二维 Ising 模型的 Monte-Carlo 方法有好几种，最常用的是单自旋翻转的 Metropolis 算法，另外还有团簇算法如 Wolff 和 Swendsen-Wang。选 Metropolis 的理由很直接：实现简单、每步只翻一个自旋、和 MATLAB 的向量化操作容易结合，而且对理解细致平衡和接受率最直观。缺点是临界点附近自相关时间急剧增大，磁化强度涨落被拉长，曲线看起来“不收敛”。Wolff 算法用团簇翻转，能把临界慢化压下去，但代码复杂度高，调试时不容易定位问题。我一般先用 Metropolis 把整条曲线跑通，确认能量和磁化强度的量级对得上，再考虑换团簇算法做临界区加密采样。如果你只是做课程设计或验证相变趋势，Metropolis 足够；如果要在 Tc 附近做高精度临界指数拟合，Wolff 更合适。

2.2 最小可跑通的 MATLAB 代码骨架

下面这段代码是我常用的骨架，格点大小 L、温度 T、蒙特卡洛步数 nsteps 都是可调参数。核心是两层循环：外层扫温度，内层做 Metropolis 翻转。为了速度，我用矩阵索引而不是逐点循环，但为了可读性，这里保留逐点写法，后面再讲向量化。

% 二维 Ising 模型 Metropolis 采样最小骨架 L = 32; % 格点边长，总自旋数 N = L*L T = 2.0; % 温度，单位取 J/kB = 1 nsteps = 20000; % 总蒙特卡洛步数 nthermal = 2000; % 热化步数，前这么多步不统计 rng(42); % 固定随机种子，保证可复现 S = 2*(rand(L) > 0.5) - 1; % 随机初始化 ±1 自旋矩阵 mag = zeros(nsteps,1); eng = zeros(nsteps,1); for step = 1:nsteps for k = 1:L*L i = randi(L); j = randi(L); % 周期边界：用 mod 处理上下左右邻居 up = S(mod(i-2,L)+1, j); down = S(mod(i,L)+1, j); left = S(i, mod(j-2,L)+1); right = S(i, mod(j,L)+1); dE = 2*S(i,j)*(up+down+left+right); % 翻转引起的能量变化 if dE <= 0 || rand() < exp(-dE/T) S(i,j) = -S(i,j); % 接受翻转 end end mag(step) = sum(S(:))/ (L*L); eng(step) = -sum(sum(S .* (circshift(S,1,1)+circshift(S,-1,1)+ ... circshift(S,1,2)+circshift(S,-1,2))))/2; end % 去掉热化段后取平均 mag_avg = mean(mag(nthermal+1:end)); eng_avg = mean(eng(nthermal+1:end)); fprintf('T=%.2f, <|M|=%.4f, <E>=%.4f\n', T, abs(mag_avg), eng_avg);

逻辑说明：dE 是翻转一个自旋后系统能量的变化，公式 dE = 2S(i,j)(邻居和) 来自 Ising 哈密顿量 H = -J Σ S_i S_j，取 J=1。接受准则用 Metropolis：dE 小于等于 0 直接接受，否则以 exp(-dE/T) 概率接受。能量统计里用 circshift 做周期边界求和，除以 2 是因为每条键被算了两次。参数说明：L 决定有限尺寸效应，L 越大临界峰越尖但计算量按 L² 增长；nsteps 要足够长，临界点附近建议至少 10⁵ 步；nthermal 一般取总步数的 10% 到 20%，具体看磁化强度是否已经稳定。rng 固定种子是为了复现，但做统计平均时应该换多组种子取平均，否则单条曲线会有偏。

2.3 把单点温度扩成整条磁化曲线

单点温度跑通后，下一步是扫温度。常见做法是 T 从 1.0 到 3.5，步长 0.1，每个温度独立跑一次。注意每个温度都要重新热化，不能把上一个温度的末态直接拿过来当初始态，除非你做的是退火模拟。下面是把上面的骨架包成函数并扫温度的写法。

function [T_list, M_list, E_list] = ising_sweep(L, T_list, nsteps, nthermal) M_list = zeros(size(T_list)); E_list = zeros(size(T_list)); for idx = 1:length(T_list) T = T_list(idx); rng(42 + idx); % 每个温度用不同种子，避免伪相关 S = 2*(rand(L) > 0.5) - 1; mag = zeros(nsteps,1); eng = zeros(nsteps,1); for step = 1:nsteps for k = 1:L*L i = randi(L); j = randi(L); up = S(mod(i-2,L)+1, j); down = S(mod(i,L)+1, j); left = S(i, mod(j-2,L)+1); right = S(i, mod(j,L)+1); dE = 2*S(i,j)*(up+down+left+right); if dE <= 0 || rand() < exp(-dE/T) S(i,j) = -S(i,j); end end mag(step) = sum(S(:))/(L*L); eng(step) = -sum(sum(S .* (circshift(S,1,1)+circshift(S,-1,1)+ ... circshift(S,1,2)+circshift(S,-1,2))))/2; end M_list(idx) = mean(abs(mag(nthermal+1:end))); E_list(idx) = mean(eng(nthermal+1:end)); end end

调用时用 T_list = 1.0:0.1:3.5，L 取 16、32、64 各跑一遍，就能看到磁化强度随温度下降的曲线，以及 L 越大曲线越陡的有限尺寸效应。这里的关键参数是 nsteps 和 nthermal 的比例，临界点附近自相关时间长，nsteps 不够会导致 M_list 在 Tc 附近跳动很大。我一般会先跑 L=16 快速看趋势，再上 L=32 或 64 做正式统计。

3. 磁化分析系统的四个核心指标：磁化强度、能量、比热和磁化率

3.1 磁化强度和能量：从原始序列到统计平均

磁化强度 M = (1/N) Σ S_i，能量 E = -Σ_{ } S_i S_j。原始序列 mag 和 eng 是每一步的值，直接取平均会受热化段影响，所以要去掉前 nthermal 步。更稳妥的做法是画一下 mag 随 step 的变化，肉眼看它什么时候稳定。如果曲线一直在漂，说明 nthermal 不够或者 nsteps 太短。统计平均时，除了均值，还要算标准差，用来估计误差棒。标准差除以 sqrt(有效独立样本数) 才是均值的标准误，但有效样本数受自相关影响，不能直接用 nsteps - nthermal。简单做法是用分块平均：把序列分成若干块，每块取平均，再看块平均的涨落。

3.2 比热和磁化率：涨落耗散定理怎么落到代码里

比热 C 和磁化率 χ 不是直接测量的量，而是从能量和磁化强度的涨落里算出来的。涨落耗散定理给出 C = (⟨E²⟩ - ⟨E⟩²) / (N T²)，χ = (⟨M²⟩ - ⟨M⟩²) / (N T)。注意这里 M 用带符号的磁化强度，不是绝对值，因为 χ 的定义涉及 M 的涨落。代码上就是在统计段里多算两个量。

% 在统计段计算比热和磁化率 E_stat = eng(nthermal+1:end); M_stat = mag(nthermal+1:end); C = (mean(E_stat.^2) - mean(E_stat)^2) / (L*L * T^2); chi = (mean(M_stat.^2) - mean(M_stat)^2) / (L*L * T); fprintf('T=%.2f, C=%.4f, chi=%.4f\n', T, C, chi);

参数说明：T 是当前温度，L*L 是总自旋数 N。比热和磁化率在 Tc 附近会出现峰值，L 越大峰越高越尖。如果你发现 C 或 chi 出现负值，通常是统计段太短或者热化不够，涨落没采样充分。另一个常见问题是 M_stat 用了 abs 之后再算 chi，这样会丢掉符号信息，导致 chi 偏小，正确做法是保留带符号的 mag 序列。

3.3 用 gnuplot 或 MATLAB 自带绘图把曲线画出来

数据跑出来后，画图是验证结果最直接的手段。MATLAB 里用 plot(T_list, M_list, ‘o-’) 就能看磁化曲线，再叠加 Onsager 解析解做对比。Onsager 给出的自发磁化强度在 T < Tc 时 m = (1 - sinh(2/T)^(-4))^(1/8)，Tc = 2/ln(1+sqrt(2)) ≈ 2.269。把解析曲线和模拟点画在一起，如果 L=32 的模拟点在 Tc 附近明显偏离解析解，那是有限尺寸效应，不是代码错。画图时建议把不同 L 的曲线放一张图，横轴 T，纵轴 |M|，能直观看到曲线随 L 增大变陡。比热和磁化率可以单独画，峰值位置随 L 增大向 Tc 靠拢。

4. 避坑与排查：Ising 模型 Monte-Carlo 模拟里最容易翻车的五件事

4.1 磁化强度在 Tc 附近剧烈跳动

现象：T 设在 2.2 到 2.4 之间，每次跑出来的 M 差别很大，曲线像噪声。原因：临界点附近自相关时间发散，Metropolis 单自旋翻转效率极低，系统在有序和无序态之间来回切换，采样不充分。解决：增加 nsteps 到 10⁵ 以上，或者换 Wolff 团簇算法；也可以做多次独立运行取平均，每次换随机种子。如果只是看趋势，把 L 减小到 16 能缓解，但会引入更强的有限尺寸效应。

4.2 能量统计出现负的比热

现象：算出来的 C 是负数，物理上不可能。原因：统计段太短，能量涨落没采样够，或者热化段没去掉，初始态的能量还在影响统计。解决：检查 eng 序列是否已经稳定，把 nthermal 调大；确认 C 公式里用的是 E_stat 的涨落，而不是把 eng 整体拿来算。另一个隐蔽原因是能量用 circshift 求和时边界重复计算，导致 E 的量级不对，涨落也跟着错。

4.3 周期边界条件写错导致磁化被钉住

现象：磁化强度一直接近 ±1，温度升高也不下降。原因：邻居索引写错，比如 mod(i,L)+1 和 mod(i-2,L)+1 搞反，导致边界自旋只和少数邻居耦合，系统被人为冻结。解决：用一个小 L（比如 4）打印邻居索引，手动验证每个格点的上下左右是否正确；或者直接用 circshift 做邻居求和，避免手写 mod 出错。周期边界是 Ising 模型的标准做法，写错会让整条曲线失去意义。

4.4 随机数种子固定导致结果不可复现的假象

现象：换了 MATLAB 版本或换了机器，同样的种子跑出来结果不一样。原因：不同 MATLAB 版本的 rand 和 randi 底层实现可能不同，randi 的序列在不同版本间不保证一致。解决：如果要严格复现，用 RandStream 显式指定生成器，比如 RandStream(‘mt19937ar’, ‘Seed’, 42)，并在代码里固定。另外，做统计平均时不要只用一个种子，至少跑 5 到 10 组不同种子取平均，否则单条曲线会有系统性偏差。

4.5 向量化写法把内存撑爆

现象：L 取 128 或 256 时，MATLAB 报内存不足。原因：一次性生成 L×L×nsteps 的三维数组，或者用 repmat 扩展邻居矩阵，内存占用按 L²×nsteps 增长。解决：不要保存每一步的完整自旋矩阵，只保存标量 mag 和 eng；邻居求和用 circshift 在每一步内做，不要预分配大矩阵。如果确实需要保存构型，每隔一定步数存一次，而不是每步都存。L=64 时，单精度存一个构型才 16KB，但存 10⁵ 步就是 1.6GB，很容易翻车。

5. 进阶技巧：用向量化和多尺度采样把临界区跑准

5.1 向量化 Metropolis：一次翻转多个自旋

逐点循环在 MATLAB 里很慢，L=32、nsteps=20000 大概要跑几分钟。向量化做法是每次随机选一批格点，批量计算 dE 和接受概率。注意批量翻转会破坏细致平衡，因为同一批里相邻格点可能互相影响。常见折中是棋盘格分解：把格点按 (i+j) 奇偶分成两组，每组内部格点不相邻，可以同时翻转而不违反细致平衡。下面是一个棋盘格向量化的核心片段。

% 棋盘格向量化：每次翻转一组不相邻格点 for step = 1:nsteps for parity = 0:1 [ii, jj] = find(mod((1:L)’ + (1:L), 2) == parity); idx = sub2ind([L L], ii, jj); up = S(mod(ii-2,L)+1 + (jj-1)*L); down = S(mod(ii,L)+1 + (jj-1)*L); left = S(ii + (mod(jj-2,L))*L); right = S(ii + (mod(jj,L))*L); dE = 2*S(idx).*(up+down+left+right); accept = (dE <= 0) | (rand(size(dE)) < exp(-dE/T)); S(idx(accept)) = -S(idx(accept)); end mag(step) = sum(S(:))/(L*L); eng(step) = -sum(sum(S .* (circshift(S,1,1)+circshift(S,-1,1)+ ... circshift(S,1,2)+circshift(S,-1,2))))/2; end

逻辑说明：parity=0 和 1 分别处理两组格点，每组内部没有最近邻关系，所以可以同时翻转。dE 是向量，accept 是逻辑向量，只翻转接受的格点。参数说明：这种写法比逐点循环快 5 到 10 倍，但要注意 find 和 sub2ind 的开销，L 很大时可以用逻辑索引代替。棋盘格分解不改变 Metropolis 的接受准则，只是把更新顺序并行化，所以细致平衡仍然满足。

5.2 多尺度采样：从高温退火到低温

临界区采样难，一个实用技巧是从高温往低温退火，每个温度用上一个温度的末态做初始态，减少热化时间。但退火会引入历史依赖，统计时要小心。更稳妥的做法是并行跑多个独立副本，每个副本从不同随机初态出发，最后取平均。MATLAB 里可以用 parfor 并行扫温度，但要注意每个 worker 的随机流独立。如果不想用并行工具箱，就手动分批次跑，把结果存成 mat 文件再合并。

5.3 用有限尺寸标定验证临界温度

判断模拟是否靠谱，一个硬指标是有限尺寸标定。对每个 L，找到比热峰值对应的温度 Tc(L)，然后按 Tc(L) = Tc(∞) + a/L 拟合，看 Tc(∞) 是否接近 2.269。如果偏离超过 0.05，说明采样或边界条件有问题。下面是一个简单的拟合片段。

% 有限尺寸标定：用比热峰值温度拟合 Tc(∞) L_list = [16 32 64]; Tc_L = zeros(size(L_list)); for k = 1:length(L_list) [~, idx] = max(C_list{k}); % C_list 是每个 L 的比热曲线 Tc_L(k) = T_list(idx); end p = polyfit(1./L_list, Tc_L, 1); Tc_inf = p(2); fprintf('拟合 Tc(∞) = %.4f，Onsager 解析值 2.2692\n', Tc_inf);

参数说明：polyfit 用 1/L 做自变量，截距就是 Tc(∞)。L_list 至少取三个点，否则拟合不可靠。如果 Tc_inf 和 2.2692 差很多，先检查比热峰值是不是被噪声干扰，可以对比热曲线做平滑再找峰。这个验证方法比单看磁化曲线更能暴露系统性问题。

5.4 我自己的习惯：先跑小尺寸，再上大尺寸，最后对解析解

我跑 Ising 模型的习惯是三步走：先用 L=8 或 16 快速跑通全流程，确认代码没有索引错误和量级错误；再用 L=32 跑正式曲线，配合多组随机种子取平均；最后把 L=32 和 L=64 的结果叠在 Onsager 解析解上，看偏离是否在有限尺寸效应预期范围内。如果 L=64 的曲线在 Tc 以下明显低于解析解，通常是热化不够或者 nsteps 太短。这套流程帮我省了很多“玄学”调试时间，也让我对 Monte-Carlo 的误差来源有更具体的感知。希望帮到你。

本文还有配套的精品资源，点击获取