简介一份面向物理、工程与统计力学学习者的 MATLAB 脚本包专门绘制麦克斯韦速度分布律的三维概率密度图像帮助读者直观理解理想气体分子速度统计规律。压缩包内共 2 个文件包含 1 个 .m 主程序脚本和 1 个 .png 示例输出图整体仅 15KB轻量易用。目前已有 437 人学习或下载适合课堂演示、课程作业、自学验证等场景具备基础 MATLAB 环境即可运行。脚本代码完整地实现了从物理参数设置、三维速度网格生成、概率密度计算到绘图的整个流程使用者可自由调整温度、分子质量等参数观察分布曲线的形态变化同时附带的示例图片为运行结果提供了直接参考便于对照检查。整份代码结构简明注释与输出配套是理解麦克斯韦分布律及 MATLAB 数据可视化的优质入门材料。1. 麦克斯韦分布律可视化为什么 MATLAB 是最顺手的那把刀麦克斯韦分布律的教材图通常是那条不对称的速率分布曲线但做气体分子模拟时真正关心的是三维速度空间中的概率密度壳层。MATLAB 在这类可视化里的优势不在公式推导而在meshgrid、slice和isosurface这些原语可以快速把高维密度函数切成可读的切片Maxwell.m这个脚本要做的就是用一步一图的方式把速度空间中的分子概率密度变成能看到形状的图形。它适合刚接触统计物理的本科生也适合需要用分布形态验证模拟结果的课题组——前者看峰位和展宽随温度怎么变后者关心归一化系数和网格分辨率是否够用。下面从理论形式、参数选择、脚本复现和排错四个角度把它拆开。2. 从公式到参数麦克斯韦分布律的数学形态与单位制选择2.1 三维速度分布与速率分布两个公式两种图像先明确概念。三维速度分布函数 f(v_x, v_y, v_z) 描述的是速度空间中单位体积的概率密度表达式是三个高斯函数的乘积f(v_x, v_y, v_z) (m / (2πkT))^(3/2) · exp( -m(v_x^2v_y^2v_z^2) / (2kT) )这个函数在速度空间的原点最大等概率面是球面沿任意方向看都是高斯形状。而速率分布 F(v) 是把速度方向积分掉之后的分布v 是速率标量F(v) 4πv^2 (m/(2πkT))^(3/2) · exp(-m v^2/(2kT))注意 F(v) 在 v0 时等于0原因是三维速度空间中半径为 v、厚度为 dv 的球壳体积是 4πv^2 dv低速时球壳体积趋近于0。很多画图失败的例子是把三维速度分布 f 的值当作 F(v) 直接画或者反过来画纵轴量级相差 4πv^2 倍自然和课本不一样。在 Maxwell.m 或自己写的脚本里首先要确定画的是哪一个如果画速度分量的截面用 f如果画速率分布曲线用 F(v)。这个选择错误是最隐蔽的坑因为两条曲线都叫麦克斯韦分布。2.2 参数表与单位制k、m、T 怎么给才不出错列一张参数速查表方便写代码时对照参数符号推荐单位常用值示例玻尔兹曼常数kJ/K1.380649e-23分子质量mkg氮气约 4.65e-26温度TK300速度分量v_xm/s网格范围 ±1200最概然速率v_pm/ssqrt(2kT/m)分子质量这里是最容易出错的。化学手册给的一般是摩尔质量比如氮气 28 g/mol要换成 kg 并除以阿伏伽德罗常数28.0e-3 / 6.02214076e23。如果直接用 28 参与计算得到的 v_p 会大好几个数量级。我用 MATLAB 写分布脚本时习惯先用下面这段代码算出 v_p再决定速度网格的上限k 1.380649e-23; % J/K m_N2 28.0e-3 / 6.02214076e23; % 氮气单分子质量kg T 300; % K v_p sqrt(2 * k * T / m_N2); % 最概然速率m/s fprintf(最概然速率 v_p %.1f m/s\n, v_p);这里的fprintf输出约 422 m/s。v_p的物理意义是 F(v) 取极大值时的速率但它同时给出了速度空间的尺度对于三维分布网格范围取 0 到 4·v_p 已经能覆盖接近 99% 的概率质量。如果温度从 300K 降到 80Kv_p 会降到 218 m/s此时网格上限还放在 1500大部分网格点落在密度极低的区域归一化和切片都会失真。所以参数选择的顺序应该是确定气体和温度 → 算 v_p → 定 vmax → 定网格点数 N。此外还要注意MATLAB 中meshgrid和ndgrid的维度含义不同。在三维情况下meshgrid生成的三个数组第一个数组随第一维行变化第二个数组随第二维列变化第三个数组随第三维页变化。如果只按教科书公式V2 Vx.^2 Vy.^2 Vz.^2计算区别不大但后续如果要用slice查看某个平面就要明白slice(Vx, Vy, Vz, f, ...)中的前三个参数必须与meshgrid输出一一对应否则切片会取错平面。这部分放到第3章展开。3. 复现 Maxwell.m网格构建、密度计算与可视化切片3.1 网格构建用 meshgrid 还是 ndgrid在三维速度空间画分布核心是构造 v_x, v_y, v_z 三个轴上的采样点再用meshgrid展开成三维数组。常见做法是N 200; % 每维点数 v linspace(-vmax, vmax, N); [Vx, Vy, Vz] meshgrid(v, v, v);这里meshgrid(v,v,v)的三个输出都是 N×N×N 数组分别对应空间中每个格点的三个速度分量。注意meshgrid输出维度顺序是 x 沿列、y 沿行、z 沿页而ndgrid相反。如果后续计算用Vx.^2Vy.^2Vz.^2两种方式没区别但用slice时meshgrid是必须的因为slice期望第一参数是 X 坐标数组且维度顺序和meshgrid一致。我一般固定用meshgrid避免在三维可视化时额外转置。网格点数 N 取多少合适200 的三维数组是 8e6 个元素每个 double 8 字节三个分量加密度共 4 个数组约 256MB普通电脑还能接受。N400 时单数组就 512MB很容易触发内存不足。工程上先取 150 验证形状再决定是否加密。3.2 计算密度注意指数计算与优化给一个可以直接放进 Maxwell.m 的密度计算段k 1.380649e-23; m 28.0e-3 / 6.02214076e23; % N2 T 300; vmax 1200; N 150; v linspace(-vmax, vmax, N); [Vx, Vy, Vz] meshgrid(v, v, v); V2 Vx.^2 Vy.^2 Vz.^2; % 逐元素平方和 f (m/(2*pi*k*T))^(3/2) .* exp(-m .* V2 ./ (2*k*T)); % 概率密度V2已经是三个方向平方之和m.*V2./(2*k*T)整体是逐元素计算。注意这里m是标量(m/(2*pi*k*T))^(3/2)也是标量和exp结果做元素乘。如果写成^而不是.^在数组上会报错或得到矩阵幂这是新手常犯的 MATLAB 错误。另外指数里温度必须用开尔文如果从摄氏温度转换要在传入前先加 273.15。这一步结束后f的数值是概率密度单位是 (m/s)^(-3)量级很小。要验证网格是否足够细可以计算数值积分dv 2 * vmax / (N - 1); integral sum(f(:)) * dv^3; fprintf(三维积分值 %.4f理想值为1\n, integral);sum(f(:))把所有格点密度相加乘以每个体积元dv^3得到积分的矩形近似。如果积分偏离 1 超过 1%通常原因是 vmax 太小截断了尾部或者 N 太小导致球壳采样不足。注意dv 2*vmax/(N-1)是因为linspace包含两端采样间隔是总范围除以 N-1。3.3 画切片与等值面slice 和 isosurface 的参数门道密度算出来后推荐先用切片看内部结构figure; slice(Vx, Vy, Vz, f, 0, [-600 0 600], 0); xlabel(v_x (m/s)); ylabel(v_y (m/s)); zlabel(v_z (m/s)); shading interp; colorbar; axis equal;这段代码在v_x0、v_y-600/0/600、v_z0的位置切出六个平面。slice的参数是(X,Y,Z,V, x_slice, y_slice, z_slice)这里第一个 0 对应 X 方向切面[-600 0 600]是 Y 方向切三个面最后一个 0 是 Z 方向切面。实际运行时会看到中心密度高、外围逐渐消失的球对称结构。shading interp消除网格线引起的色块边界axis equal保证三个方向单位长度一致否则球壳会被拉伸成椭球。如果觉得切片不够直观用等值面把指定密度阈值的曲面拉出来threshold max(f(:)) * 0.2; % 取峰值的20%作为等值面阈值 isosurface(Vx, Vy, Vz, f, threshold);max(f(:))*0.2这个阈值是相对值因为 f 的绝对量级随 m 和 T 变化很大绝对阈值难估。网络搜索里常有人问为什么isosurface画出来是空图原因基本都是阈值给得比实际最大值还大。先用max看一下峰值再定阈值比凭感觉给一个 1e-6 之类的数可靠得多。网格点数直接决定内存和运行耗时我习惯参考下表选择 NN四个数组总内存运行耗时三维积分精度50约 4 MB极快1.02100约 32 MB快1.005150约 108 MB中等1.001200约 256 MB较慢1.00054. 参数扫描温度、分子量与分布曲线形态的关系4.1 从最概然速率看懂参数影响分布曲线有两个关键特征峰位和展宽。峰位由 v_p sqrt(2kT/m) 决定温度越高或分子质量越小峰值对应速率越大。展宽则由整个高斯指数项决定温度高时分布更平缓。给一个速查表气体分子质量(kg)300K v_p (m/s)600K v_p (m/s)1000K v_p (m/s)He6.65e-27111815812040N24.65e-26422597771CO27.31e-26337477615这个表能看出温度翻4倍v_p 翻1倍因为 v_p 与 T 的平方根成正比。做参数扫描时也按这个规律选温度点比如 100、200、400、800K等比分布比等差分布更能体现曲线形态变化。4.2 一维速率分布的多温度对比代码这里给一个能直接跑出对比图的脚本段m 4.65e-26; % N2 分子质量 k 1.380649e-23; v 0:5:1600; % 速率从0开始禁止负值 T_list [100, 300, 600]; for i 1:length(T_list) T T_list(i); F 4*pi*v.^2 .* (m/(2*pi*k*T))^(3/2) .* exp(-m*v.^2/(2*k*T)); plot(v, F/1e-6, LineWidth, 1.5); hold on; end xlabel(速率 v (m/s)); ylabel(F(v) (×10^{-6} s/m)); legend(100K,300K,600K);代码里F/1e-6只是缩放纵轴避免数值太小在图形上显示为 10^-6 量级实际数据如果要用保留原始F即可。注意这里用v.^2元素运算4*pi*v.^2是每一项乘以对应速率平方。为什么速率从0开始因为 F(v) 的定义域是 v0负速率没有物理意义。如果习惯性写v linspace(-1000,1000)画出来的图形会多出左侧对称的一半不是麦克斯韦分布。4.3 归一化、坐标范围和内存的坑第一归一化混淆。三维分布 f 和速率分布 F(v) 的归一化条件不同f 对 v_x v_y v_z 三重积分等于1F(v) 对 v 从0到∞积分等于1。如果在三维脚本里用sum(f(:))直接当作粒子数会得到 1/dv^3 量级的大数必须乘上体积元。同理一维速率分布画出来后用trapz(v, F)验证积分是否接近1不要用sum(F)。第二坐标范围。温度 300K 的氮气v_p422速率分布的主要区间在 0~1200。如果把v设成 0:5:300只能看到上升沿之前的极小部分曲线看起来像指数衰减。反过来网格上限太大而 N 不够的时候每个格点之间的间隔变大等值面会变成锯齿状。我的经验是vmax 4*v_pN 取 150~200这样分布主体和尾部都有足够分辨率。第三内存。三维meshgrid对 N 的内存增长是立方。N200 时四个 double 数组约 256MB加上图形对象还能接受N400 时约 2GB在旧电脑上直接卡死。想要更高分辨率可以用single精度存储Vx,Vy,Vz或者只计算V2而不是三个分量数组。如果内存仍不够就降为二维切片固定一个速度分量为0用contourf画二维截面代码在 3.3 节基础上把slice换成contourf即可。5. 进阶验证用积分和动画检验 Maxwell.m 的归一化精度5.1 把积分验证做成一个可复用函数开发Maxwell.m时我习惯在脚本开头放一个自检函数避免参数改坏后曲线形态变了却不知道是网格问题还是公式问题。下面这个函数接收网格点数、温度和气体名称输出数值积分值和推荐网格范围function [integral, vmax_recommend] check_integral(N, T, gas) k 1.380649e-23; switch gas case N2, m 28.0e-3 / 6.02214076e23; case He, m 4.0026e-3 / 6.02214076e23; case CO2, m 44.01e-3 / 6.02214076e23; end v_p sqrt(2*k*T/m); vmax_recommend 4*v_p; v linspace(-vmax_recommend, vmax_recommend, N); [Vx, Vy, Vz] meshgrid(v, v, v); V2 Vx.^2 Vy.^2 Vz.^2; f (m/(2*pi*k*T))^(3/2) * exp(-m*V2/(2*k*T)); dv 2*vmax_recommend/(N-1); integral sum(f(:)) * dv^3; fprintf(N%d, T%gK, vmax%.1f, integral%.4f\n, N, T, vmax_recommend, integral); end这个函数把气体摩尔质量到单分子质量的转换封装在switch中避免每个脚本重复写28e-3/6.022e23。运行时用check_integral(150, 300, N2)输出积分值一般在 1.001 左右。如果 N 加到 250积分会更接近 1但耗时翻倍。这个验证能发现两类问题一是vmax不足导致积分小于 1二是温度极低时分布尾部在默认网格下采样不足。5.2 把温度扫描做成 GIF 动画静态图只能看一条曲线动画能更清楚展示峰位和峰高的联动关系。简单写法如下v 0:5:1600; frames cell(1, 11); for T 100:50:600 F 4*pi*v.^2 .* (m/(2*pi*k*T))^(3/2) .* exp(-m*v.^2/(2*k*T)); plot(v, F/1e-6, LineWidth, 1.5); ylim([0, 3]); xlim([0, 1600]); title(sprintf(T %d K, T)); frames{end1} getframe(gcf); endframes需要预先分配或动态扩展温度从 100 到 600 每次加 50共 11 帧。sprintf把温度写进标题方便逐帧查看。导出 GIF 时新版 MATLAB 可以直接用exportgraphics(gcf,maxwell.gif,Resolution,100)或者用imwrite循环写帧。如果想把温度序列改成更密集的100:10:200动画前半段的变化会比后半段明显得多这是 v_p 与 T 平方根关系最直观的演示。本文还有配套的精品资源点击获取