★ FEATURED ARTICLE 简介这份资源面向凝聚态物理、材料科学方向的学习者与科研人员提供一套用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)*(updownleftright); % 翻转引起的能量变化 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(nthermal1:end)); eng_avg mean(eng(nthermal1: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取 J1。接受准则用 MetropolisdE 小于等于 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)*(updownleftright); 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(nthermal1:end))); E_list(idx) mean(eng(nthermal1:end)); end end调用时用 T_list 1.0:0.1:3.5L 取 16、32、64 各跑一遍就能看到磁化强度随温度下降的曲线以及 L 越大曲线越陡的有限尺寸效应。这里的关键参数是 nsteps 和 nthermal 的比例临界点附近自相关时间长nsteps 不够会导致 M_list 在 Tc 附近跳动很大。我一般会先跑 L16 快速看趋势再上 L32 或 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(nthermal1:end); M_stat mag(nthermal1: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(1sqrt(2)) ≈ 2.269。把解析曲线和模拟点画在一起如果 L32 的模拟点在 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 在每一步内做不要预分配大矩阵。如果确实需要保存构型每隔一定步数存一次而不是每步都存。L64 时单精度存一个构型才 16KB但存 10⁵ 步就是 1.6GB很容易翻车。5. 进阶技巧用向量化和多尺度采样把临界区跑准5.1 向量化 Metropolis一次翻转多个自旋逐点循环在 MATLAB 里很慢L32、nsteps20000 大概要跑几分钟。向量化做法是每次随机选一批格点批量计算 dE 和接受概率。注意批量翻转会破坏细致平衡因为同一批里相邻格点可能互相影响。常见折中是棋盘格分解把格点按 (ij) 奇偶分成两组每组内部格点不相邻可以同时翻转而不违反细致平衡。下面是一个棋盘格向量化的核心片段。% 棋盘格向量化每次翻转一组不相邻格点 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).*(updownleftright); 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逻辑说明parity0 和 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(∞) %.4fOnsager 解析值 2.2692\n, Tc_inf);参数说明polyfit 用 1/L 做自变量截距就是 Tc(∞)。L_list 至少取三个点否则拟合不可靠。如果 Tc_inf 和 2.2692 差很多先检查比热峰值是不是被噪声干扰可以对比热曲线做平滑再找峰。这个验证方法比单看磁化曲线更能暴露系统性问题。5.4 我自己的习惯先跑小尺寸再上大尺寸最后对解析解我跑 Ising 模型的习惯是三步走先用 L8 或 16 快速跑通全流程确认代码没有索引错误和量级错误再用 L32 跑正式曲线配合多组随机种子取平均最后把 L32 和 L64 的结果叠在 Onsager 解析解上看偏离是否在有限尺寸效应预期范围内。如果 L64 的曲线在 Tc 以下明显低于解析解通常是热化不够或者 nsteps 太短。这套流程帮我省了很多“玄学”调试时间也让我对 Monte-Carlo 的误差来源有更具体的感知。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
