二维Ising模型蒙特卡洛磁化分析：Metropolis算法与MATLAB实现

简介：基于Monte-Carlo模拟的二维Ising模型磁化分析系统，是一份使用MATLAB实现的数值模拟工具，面向凝聚态物理、材料科学等领域的研究人员和学生。它可在不开展复杂物理实验的情况下，预测铁磁材料在不同温度下的磁性能，并帮助用户理解磁化强度、能量等热力学量随温度的演化及相变机制。资源包共2个文件，包含MATLAB主程序脚本与Markdown说明文档，压缩包仅4KB，体量轻、结构清晰，方便快速部署与阅读。主程序基于Metropolis算法完成Monte-Carlo迭代：随机选取自旋、计算翻转前后能量变化，按能量判据接受或拒绝翻转，最终逼近热力学平衡；可设置网格尺寸、初始温度等参数，输出平均磁化强度与能量随温度变化曲线，并展示不同温度下自旋网格的排布可视化，直观呈现临界温度附近的磁序突变。配套说明文档梳理了模型原理、算法流程、参数设置与结果解读，便于二次开发和教学演示；模拟结果与铁磁材料的实际物理现象相符，为研究相变规律提供了便捷参考。目前已有105人学习浏览，适合统计物理课程、数值模拟入门及磁相变研究的辅助工具。

1. 二维Ising模型Monte-Carlo磁化分析：第一版就能复现Tc≈2.269

这套基于MATLAB的二维Ising模型Monte-Carlo模拟磁化分析系统，核心就一件事：用Metropolis算法在L×L的格点上做自旋翻转采样，再把统计平均的温度依赖关系画成磁化强度-温度曲线。我第一次跑通时最直接的感受是，磁化强度曲线在临界温度附近肉眼可见地跌落，不用做任何拟合就能看出相变发生，这种“计算复现物理”的直观冲击比看任何教材公式都强。这套系统适合两类人：一类是统计物理、计算物理方向做课程设计或毕业设计的学生，另一类是刚接触蒙托卡洛模拟、想在MATLAB里跑通一个完整体系的工程师。它最大的价值不是算法多难，而是把“随机初态、热平衡丢弃、MCS采样、涨落统计”这一整条链路揉成了一组能直接改参数的函数文件。

2. Metropolis核心逻辑：哈密顿量、能量差与翻转概率怎么落地

2.1 二维Ising模型的哈密顿量与单自旋翻转的能量差

二维Ising模型里，每个格点上的自旋s只取+1或−1，体系哈密顿量写成

H = −J·Σ_⟨i,j⟩ s_i·s_j − h·Σ_i s_i

其中⟨i,j⟩表示最近邻自旋对，二维正方格点上每个自旋有上下左右四个近邻；J是耦合常数，J>0时体系倾向于让相邻自旋同向排列，这正是铁磁序的来源；h是外磁场。模拟里我只保留最近邻相互作用，不做次近邻扩展，这样能和Onsager的精确解直接对照。

Metropolis更新时需要对单个自旋翻转前后的能量差做精确计算，公式很短：

ΔE = 2·J·s·(邻居自旋和) + 2·h·s

这个公式是从“翻转前后哈密顿量相减”推出来的：翻转一个自旋，只有包含它的那些项变号。我在第一版代码里图省事，直接对整个格点重算总能量再做差，L=16时还能忍，L=64后同一个温度扫描慢得像龟爬。改成单自旋局部能量差后，每个MCS的复杂度从O(L²)量级降到一个常数级别。实际计算邻居和时，周期性边界是最容易写错的地方，我习惯的写法是：

% 周期性边界下的四个近邻和，L为格点尺寸 nb = lattice(mod(x-2, L)+1, y) ... % 上方邻居 + lattice(mod(x, L)+1, y) ... % 下方邻居 + lattice(x, mod(y-2, L)+1) ... % 左方邻居 + lattice(x, mod(y, L)+1); % 右方邻居

这里的mod索引写法值得单独说：MATLAB的mod(0,L)=0，所以取“绕回边界”时要加1，mod(x-2,L)+1在x=1时指向第L行，在x=2时指向第1行，逻辑正好对得上。我见过不少翻车案例是把这里写成if边界判断，结果四个角的邻居数变成3个，相当于默认开了自由边界，磁化曲线整体向低温方向偏移。验证边界对不对有个笨办法：统计每个格点的邻居数，必须是4，不能出现3或2。另外注意，局部能量差公式里外场项也要乘2，很多人漏掉h那项，在h≠0的模拟里磁滞回线会错得离谱。

2.2 Metropolis接受概率与“只能下降”的错误直觉

Metropolis算法的核心是一个接受概率：

p_accept = min(1, exp(−β·ΔE))

其中β=1/(k_B·T)。温度T的单位在模拟里按k_B=1处理，所以T是J/k_B的无量纲倍数，临界温度T_c≈2.269就是这个约定下的值。ΔE≤0的翻转直接接受，系统能量只会下降；ΔE>0时以概率exp(−βΔE)接受，允许短暂的能量上升，让体系能越过能垒，从局部极小里逃出来。

这里有个新手特别容易犯的直觉错误：认为既然是找基态，就应该拒绝一切能量上升的翻转。如果真这么干，低温模拟会直接掉进完全有序态，但高温模拟也回不到无序态，因为体系永远找不到从有序到无序的路径。蒙特卡洛模拟要的是“按照玻尔兹曼分布采样”，不是找能量极小值。在实际执行时，我会把接受判断写成：

if dE <= 0 || rand < exp(-beta * dE) lattice(x, y) = -s; % 接受翻转 end

注意逻辑顺序：dE>0时exp(−β·dE)必然小于1且大于0，当随机数rand小到能落在这个概率内就接受；dE很大时exp可能下溢成0，翻转被拒绝，这在物理上也是正确的。MATLAB里短路的||运算符能避免不必要的exp计算，这点比直接写rand < exp(...)再单独判断要快一点，温度越低越明显。

2.3 MCS步数约定、burn-in与物理量的涨落统计

一个MCS（Monte-Carlo Step）的定义是：对格点上所有自旋平均尝试一次翻转，也就是L²次Metropolis尝试。我实现的版本里直接对每个格点做一次扫描，这样一遍循环就是一个MCS，统计步数时不会产生歧义。体系从随机初态出发，前几个MCS能量会剧烈下降，这是热化过程，这段数据不能参与统计，必须丢掉的步数叫burn-in步数。

我在参数里默认给到1000步热平衡丢弃，然后继续采样5000步。热平衡判断还有一个土办法：输出能量随时间变化的曲线，如果尾部稳定在一个水平上下波动，没有单调漂移，就说明平衡了。统计量方面，磁化强度m是每自旋平均的 Σs_i/L²；磁化率和比热直接由涨落给出，不用做有限差分，这一点留到第3章代码里具体展开。模拟实践里，低温下从随机初态出发很容易形成畴结构导致m偏低，所以我在主脚本里所有初态都用rand(L)>0.5来随机生成，而不用全+1的冷启动，除非刻意要模拟单畴。

3. MATLAB代码实现：初始化、MCS更新与磁化统计的完整流程

3.1 主脚本与函数文件的组织方式

这套系统按“主脚本调三个子函数”的方式组织，主脚本负责温度扫描、统计汇总和绘图，三个子函数分别承担初始化、一个MCS的更新、以及单次构型的测量。拆成这样主要是为了排查方便：出错时能直接定位到某一个函数，而不是在一坨脚本里翻来翻去。主脚本的框架长这样：

%% ising_main.m 主脚本 clear; clc; L = 32; % 格点尺寸（正方形） T_list = 0.5:0.05:5.0; % 温度扫描范围，单位 J/k_B J = 1; % 近邻耦合强度 h = 0; % 外磁场，本次扫描固定为0 n_burn = 1000; % 热平衡丢弃步数（MCS） n_mcs = 5000; % 采样步数（MCS） % 预分配输出数组 M_avg = zeros(size(T_list)); Sus = zeros(size(T_list)); Cv = zeros(size(T_list)); for tidx = 1:numel(T_list) beta = 1 / T_list(tidx); lattice = ising_init(L); % 随机初态 for step = 1:n_burn lattice = ising_mcs(lattice, beta, J, h, L); end sum_m = 0; sum_m2 = 0; sum_e = 0; sum_e2 = 0; for step = 1:n_mcs lattice = ising_mcs(lattice, beta, J, h, L); [m, e] = ising_measure(lattice, J, h); sum_m = sum_m + m; sum_m2 = sum_m2 + m^2; sum_e = sum_e + e; sum_e2 = sum_e2 + e^2; end M_avg(tidx) = abs(sum_m / n_mcs); % 对|M|取平均 Sus(tidx) = beta * L*L * (sum_m2/n_mcs - (sum_m/n_mcs)^2); Cv(tidx) = beta^2 * L*L * (sum_e2/n_mcs - (sum_e/n_mcs)^2); end

这段代码里有几个参数是刻意留出来供调试验证的：n_burn和n_mcs直接影响统计涨落大小，L决定有限尺寸效应的强度。我通常先跑L=16的快速档把Tc的大致范围摸出来，再在临界区附近用更大的L和更细的温度步进加密。统计公式这边的血泪经验是：m是每自旋的磁化，磁化率χ=β·N·var(m)，这里N=L²必须乘上去；比热C=β²·N·var(e)同理，e是每自旋能量。如果把N漏乘，峰位虽然不变但数值完全对不上解析解的幅度，会让人误以为代码写错了。另外，对M取|M|再平均是有限尺寸模拟的常规操作，因为有限体系磁化在正负间来回翻转，直接平均会得到接近0的假象，只有|M|才能反映有序度。

3.2 初始化函数与随机数映射

初始化函数很短，但有一个常见误区是使用randn高斯随机数给自旋赋值。Ising模型自旋是离散的±1，用高斯分布还得再做阈值化，纯属多余。正确做法是均匀分布随机数与0.5比较，再映射到±1：

% ising_init.m function lattice = ising_init(L) % L: 格点行数和列数，返回 L×L 的 ±1 矩阵 lattice = rand(L) > 0.5; % 每个格点独立且以50%概率为真 lattice = lattice * 2 - 1; % true -> +1, false -> -1 end

这里rand(L)生成L×L的[0,1)均匀分布随机数，>0.5产生布尔矩阵，乘2减1后映射到+1/−1。为什么用随机初态而不是全+1？在高T区域，从全+1出发系统会表现出虚假的短时间有序，磁化衰减曲线的前段会带记忆效应。随机初态在高T下能更快进入真正的平衡。T很低时随机初态反而容易形成畴，所以我在主脚本里保留了随机初态，但建议在需要精确低温结果时改成热启动，这一点第5章排查里会再讲。

3.3 一个MCS的Metropolis更新函数

更新函数是本系统的核心计算单元，采用逐个格点扫描的方式完成一个完整MCS。函数定义如下：

% ising_mcs.m function lattice = ising_mcs(lattice, beta, J, h, L) % 对格点顺序扫描一遍，完成一个MCS % lattice: L×L的±1矩阵 % beta: 1/(k_B*T)，热力学 beta % J, h: 耦合常数与外场 for x = 1:L for y = 1:L s = lattice(x, y); % 周期边界下四个近邻自旋的数值和 nb = lattice(mod(x-2, L)+1, y) ... + lattice(mod(x, L)+1, y) ... + lattice(x, mod(y-2, L)+1) ... + lattice(x, mod(y, L)+1); dE = 2 * J * s * nb + 2 * h * s; if dE <= 0 || rand < exp(-beta * dE) lattice(x, y) = -s; % 翻转自旋 end end end end

这段代码的要点有三个。第一，nb的计算用mod移位实现周期性，不在边界处加if判断，保证四个近邻永远齐全。第二，dE的算法直接对应第2章推导的局部能量差，没有重新扫描全晶格，速度差别在大L下是数量级级别的。第三，顺序更新意味着当前格点翻转后立刻影响后面格点的邻居和，这与同步更新略有差别，但对二维Ising模型的平衡态结果没有影响，收敛速度反而更快。如果要改成随机选取L²个格点的“随机序更新”，只需把双循环改成for step=1:L*L再用randi(L)取坐标，二者统计结果一致。需要注意MATLAB里exp和rand都是向量化友好的函数，但在双循环里每次调用开销不小，如果追求性能可以在临界区附近减少采样频率，而不是优化这个循环内部。

3.4 测量函数与能量统计的细节

测量函数返回每自旋的磁化强度m和每自旋能量e。能量计算的坑点在于最近邻相互作用每条键被两个自旋分别计入一次，求和后必须除以2，否则能量曲线会系统性偏高，比热峰的形状也会偏胖：

% ising_measure.m function [m, e] = ising_measure(lattice, J, h) % 返回每自旋平均磁化强度m和每自旋平均能量e % 温度相关量通过外部beta的涨落统计获得 L = size(lattice, 1); m = mean(lattice(:)); % 每自旋磁化 nb = zeros(L, L); for x = 1:L for y = 1:L nb(x, y) = lattice(mod(x-2, L)+1, y) ... + lattice(mod(x, L)+1, y) ... + lattice(x, mod(y-2, L)+1) ... + lattice(x, mod(y, L)+1); end end e = -J * sum(sum(lattice .* nb)) / 2 - h * sum(lattice(:)); e = e / (L * L); end

这个函数每次测量都要重新算一遍邻居和，采样步数多的时候会占用不少时间。我一般不会在耗时的同时去优化它，因为5000步采样已经能得到平滑的曲线，再加大步数收益不大。如果确实追求速度，可以在更新函数返回时顺便把邻居矩阵带出来，但那样代码耦合度变高，不利于读者拆开复现。主脚本里的磁化率和比热完全靠var(m)与var(e)计算，这也是Metropolis-plus-涨落统计的标准做法，比用磁化曲线数值微分求χ要稳定得多，尤其是临界区附近数据有涨落，数值微分会产生明显噪声。

4. 温度扫描实验：三档参数配置与相变峰定位

4.1 参数配置表：不同目标下的推荐组合

温度扫描实验的参数配比直接影响模拟可信度和耗时。我按使用场景整理了三档配置，先跑快速档确认物理图像，再决定是否上精细档。这里给出我在复现时用的默认参数组合：

快速档在普通桌面机上单温度运行约1~2秒，整个扫描约1分钟；标准档单温度10秒左右，50个温度点约8分钟；精细档耗时按照L²增长，单温度会到1分钟以上，临界区附近因为自关联变长还需要额外加采样。我强烈建议先跑快速档看全局趋势，确认Tc落在预期位置附近后再把温度轴改到1.8到2.8之间做精细扫描，这样能省掉大量无意义的高低温数据。另一个习惯是温度步进在临界区缩小到0.02，因为χ峰的形状对步进很敏感，步进太大会直接错过峰值位置。

4.2 自旋组态可视化：直接从格点图看出有序与无序

温度扫描过程中，把特定温度下的自旋组态用imagesc画出来，比任何数理统计都直观。在标准档里我会在T=2.0、T=2.269、T=3.0三个温度点保存组态图像：

% 自旋组态可视化片段 figure; subplot(1,3,1); imagesc(lattice_T200); axis square; colormap(gray); title('T = 2.00，铁磁畴结构'); subplot(1,3,2); imagesc(lattice_Tc); axis square; colormap(gray); title('T ≈ 2.269，临界涨落'); subplot(1,3,3); imagesc(lattice_T300); axis square; colormap(gray); title('T = 3.00，顺磁无序');

灰度图里+1格点显示为白色，−1格点为黑色。T=2.0时能看到大片同色区域连成块，这是铁磁畴；T≈2.269时黑白斑块尺度变大并且形态支离破碎，这正是临界涨落导致的关联长度发散；T=3.0时黑白格点均匀混合，没有明显大块，说明体系进入顺磁相。有一个肉眼可直接判断系统是否正常的技巧：T远低于Tc时如果看到蛛网一样的细长畴壁，说明从随机初态出发后体系被冻结在亚稳态，需要更长的burn-in或者换热启动。

4.3 用χ峰与比热峰精确定位相变临界点

|M|-T曲线用来判断相变存在很直观，但要报告一个数值化的Tc，我通常直接去找磁化率χ的峰值位置。χ峰比|M|曲线拐点更尖锐，对临界涨落响应更明显，在小体系里定位也更稳定。定位代码非常简单：

% 从磁化率序列定位峰值温度 [~, idx] = max(Sus); Tc_est = T_list(idx); fprintf('峰值温度估计: %.3f\n', Tc_est);

需要注意，L=16的χ峰位置通常会落在2.5附近而不是2.269，这是因为有限体系把临界点往高温方向推了；L=32时峰位落在2.4左右，L=64时接近2.35。峰位随L增大逐渐逼近2.269，这就是有限尺寸标度的直观表现，具体外推方法放在第6章。如果你跑出来的L=32峰位在2.0以下，优先怀疑边界条件错误，检查办法是把格点角落的邻居数打印出来，四个邻居一个都不能少。临界区还有一个特点：χ对采样步数极其敏感，临界点上自关联时间发散，5000个MCS的独立样本数锐减，峰高会明显偏低。解决方法是临界区单独加采样到10000步甚至20000步，或者对同一温度做多次独立运行再平均。

5. 避坑与排查：磁化曲线不收敛的五类高发原因

5.1 三类常见异常现象与根因

我在复现和调试这套模拟系统的过程中，整理了五类最常踩的坑，按“现象→原因→解决”的格式记录下来。

现象一：高温区|M|平均不归零，稳定在一个0.2左右的小平台，温度越高也不掉下去。

原因：高温下体系磁化确实围绕0涨落，但取绝对值后再平均会产生正偏差；更常见的是采样步数太少，体系没有充分遍历磁场反转的两种符号。小L体系里磁化会在正负之间整段翻转，如果你恰好截断在某一符号的区间，|M|平均值就是正偏差。

解决：先加长采样到10000个MCS；如果平台仍不消失，对每个温度跑5次独立实验，把每次末段磁化按符号分垛，以多数符号为准做“有序侧”平均。我在标准档里采用abs处理主序参量，但心里清楚这只是权宜之计，报告“平均磁化”时反而应该说明用的是|M|还是符号对齐后的M。

现象二：磁化率χ峰位置每次运行都不一样，偏差超过0.1，临界区尤其严重。

原因：临界点附近自关联时间τ发散，有限采样长度下χ涨落被严重低估，同时不同初态轨迹确实会带来统计差异。这不是代码随机数出问题，而是物理系统本身的临界慢化。

解决：固定rng种子后做同种子对比，把“实现代码是否稳定”和“统计是否收敛”分开判断。代码里每次跑之前加一行rng(2024)，把种子写死；确认实现稳定后，再对临界区温度点单独把采样加到20000步。临界区的物理涨落是真实的，多跑几组取平均还能给χ峰高度画误差棒。

现象三：计算得到的Tc明显偏低，L=16时只有2.0左右，L=32也只有2.2。

原因：周期性边界条件写错。如果格点角落的邻居数少于4，相当于体系实际是开放边界，边界上的键缺失会让体系能量偏低，相变更容易发生，Tc自然左移。

解决：写一段临时脚本统计每个格点的邻居数，逐一打印角落值，必须是4。另一种隐蔽错误是邻居索引把mod用成了rem，MATLAB的rem负数结果与mod不同，也会造成同样的边界割裂。养成检查mod(-1,L)结果的习惯即可。

现象四：低温区磁化强度卡在0.9左右上不去，组态图里出现明显的反铁磁条纹或畴壁。

原因：从随机初态出发，低温下自旋翻转被强烈抑制，体系停留在大块磁畴的亚稳态里，正向磁化区与反向磁化区相互抵消，整体|M|上不去。

解决：低温段改成热启动，把上一个温度T+ΔT收敛的终态作为当前温度的初态，这相当于真实实验里的缓慢降温过程。在温度扫描循环里，把上一温度的lattice直接传给下一温度，只在最高温处用随机初态。这个方法对临界区以外区域非常有效。

现象五：下载的脚本打开后中文注释是乱码，报错信息里的中文字符也变成菱形问号。

原因：较老版本的MATLAB默认用GBK编码读取.m文件，而脚本文件用UTF-8保存，中文注释全被误读。

解决：在MATLAB的预设项里设置编辑器编码为UTF-8；或者在脚本文件上右键另存为，选择UTF-8编码。这个问题和模拟代码逻辑无关，但会让首次复现的体验极差，建议拿到脚本第一件事就是统一编码。

5.2 排查顺序：先输出能量曲线，再看磁化曲线

遇到任何“模拟结果不对”的反馈，我的排查顺序永远是先看单温度能量曲线，再看整条磁化曲线。两者结合能快速区分“没热平衡”和“逻辑写错”。检查能量曲线的代码片段：

% 单温度能量收敛检查 T_test = 2.27; % 临界区内 beta = 1 / T_test; lattice = ising_init(32); E_rec = zeros(1, 3000); for k = 1:3000 lattice = ising_mcs(lattice, beta, 1, 0, 32); [~, e] = ising_measure(lattice, 1, 0); E_rec(k) = e; end plot(E_rec); xlabel('MCS'); ylabel('每自旋能量');

如果曲线前几百步有明显单调下降，之后稳定在一个水平线附近波动，说明热平衡基本达成；如果曲线到最后还在持续下降，说明burn-in不够，需要把丢弃步数往上加；如果曲线呈现周期性的锯齿形状，那多半是代码里意外地在每个MCS内重置了随机数种子，导致“模拟”实际是一系列重复轨迹的拼接。能量曲线尾部波动的幅度还能顺势估计自关联时间，波动衰减得很慢说明临界慢化明显，需要把采样步数拉长或者干脆用多链并行。做完这步再看磁化曲线，如果能量正常但磁化失稳，问题几乎都出在初态选择和统计口径上，而不是算法框架。

6. 进阶技巧：有限尺寸标度、固定随机种子与解析解对拍

6.1 用有限尺寸标度外推真实Tc

小格点的χ峰位置会系统性偏离真临界点，把不同L的峰值位置外推到L→∞是标准的校正方案。二维Ising模型的临界指数ν=1，因此峰值位置满足T_c(L)=T_c(∞)+a/L。实际操作时，对不同L分别跑温度扫描，记录各自的峰值温度，然后做1/L的线性拟合，截距就是外推的T_c(∞)：

L_set = [16 24 32 48 64]; Tc_L = zeros(size(L_set)); for k = 1:numel(L_set) % 假设已经通过温度扫描得到 T_grid 与 Sus [T_grid, Sus] = run_scan(L_set(k)); [~, idx] = max(Sus); Tc_L(k) = T_grid(idx); end p = polyfit(1 ./ L_set, Tc_L, 1); Tc_inf = p(2); % 截距即外推临界温度

这里每个L的跑法应该用临界区精细档，温度步进至少0.02，否则峰位定位误差会盖过外推本身。拟合完成后与Onsager精确解2.269对比，误差在2%以内说明实现基本正确。L取偶数比奇数好，因为奇数格点无法保证正反磁化态的对称性，统计涨落会不对称。

6.2 固定随机种子、初态约定与并行加速的边界

可复现性是这类模拟的命根子。我在主脚本开头固定写rng(2024)，这样每次运行得到完全相同的温度和磁化曲线，方便与别人对照。排查代码逻辑错误时，固定种子还意味着同一次运行的不同阶段可以直接差分比较。在临界区附近固定种子意义更大，因为临界慢化会导致结果对初态和噪声极其敏感。

关于并行加速，常见做法是给温度扫描开MATLAB parallel pool，但临界区自关联时间很长，单个温度的模拟本身就需要长轨迹，并行能节省的是“跑多组初值取平均”的时间，而不是单链的时间。如果你的环境中提示MATLAB no parallel pool不可用，直接用串行循环跑，L=32的标准档8分钟完全可接受，不必为了省这点时间引入并行复杂度。

6.3 与解析解对拍，以及我现在的固定流程

对拍验证是判断整套系统是否可信的最后一关。二维Ising模型提供了罕见的精确解：T_c=2.269J/k_B，零场下磁化强度按临界指数β=1/8趋近于0，比热在临界处呈现对数发散。我会把模拟的χ峰外推值与精确解比较，同时确认M-T曲线在低温端形状接近β=1/8的幂律趋势即可。还要检查能量曲线：T=0.5时每自旋能量应逼近−2J，因为完全有序态每条键贡献−J/2，每个格点4条键除2后正好−2J。如果这个数值对不上，优先检查能量求和是否漏了除以2。

从那以后，我每次跑温度扫描都强制自己先走一遍固定流程：固定随机种子→跑快速档确认峰值区间→在临界区加密温度步进→用能量曲线确认热平衡→保存|M|与χ的原始序列→最后做不同L的外推对拍。这套下来基本杜绝了“曲线看起来合理但实际是错的”这类玄学问题。这套MATLAB实现里的函数文件都不长，拆开重拼也很容易，关键是每一步的统计口径和边界处理都要自己重新确认一遍。希望帮到你。

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