简介针对普朗克定律中不同温度下黑体光谱辐射强度随波长分布的可视化需求Matlab代码包提供了可直接运行的绘图方案。面向热辐射、光学或物理相关课程的本科生与研究生也适合需要快速绘制多曲线对比图的科研人员。压缩包内含3个文件包括2个.m源码文件和1个png结果图源码实现普朗克公式计算并支持在同一坐标系中绘制多条温度曲线且采用不均匀刻度图展示波长的宽动态范围便于观察峰值位移规律。结果图可直接作为论文或报告的插图参考。包体总大小约227KB内容精简无多余依赖。当前已有8371人学习下载说明该工具在相关教学与科研场景中具有较高的实用价值。通过运行示例读者可以掌握Matlab实现物理公式曲线绘制的基本流程并灵活调整温度参数生成自定义光谱分布图。1. 普朗克定律为什么温度越高曲线峰值越往短波跑黑体辐射的横轴是波长纵轴是光谱辐射强度。学过普朗克定律的人都知道结论——温度升高时整条曲线抬高峰值位置往短波方向移动。可一旦亲手在 MATLAB 里把这簇曲线画出来大概率会翻车300K 的曲线贴着零轴6000K 的曲线却冲上天线性坐标下根本没法同时看清。这里的核心障碍不是公式复杂而是“峰值移动”和“总量四次方增长”是两个独立机制前者来自指数项里的温度后者来自 λ⁵ 分母。下面从公式的来历、数值实现的细节、多曲线图与不均匀刻度画法一路拆到那些最容易踩坑的地方读完你会真正理解为什么黑体辐射曲线长这样。2. 普朗克公式的来历一条内插曲线如何终结紫外灾变2.1 紫外灾变经典物理在短波端的裂缝19 世纪末实验测出的黑体辐射曲线长得很“规矩”长波端强度随波长增大而下降短波端快速趋近于零中间有一个清晰的峰。但经典统计力学给出的瑞利-金斯公式却在长波端和实验吻合、在短波端彻底发散——波长越短辐射强度按 λ⁻⁴ 无限上升。这个矛盾被物理学家称为“紫外灾变”它暴露了“能量连续”这个假设的缺陷。普朗克当时的思路不是推翻经典理论而是做数学上的“内插”找到一个解析式让它在长波极限退化成瑞利-金斯在短波极限退化成当时已知的维恩近似中间平滑衔接。这个式子就是今天看到的普朗克公式通常写成光谱辐出度形式M_λ(λ, T) (2πhc² / λ⁵) · 1 / (e^{hc/(λkT)} − 1)这里 h 是普朗克常数c 是光速k 是玻尔兹曼常数。理解这个公式的关键在于分母里的指数项h c / (λ k T) 是一个无量纲数温度越高、波长越长它就越小当它远小于 1 时e^x − 1 ≈ x公式退化为 2πc k T / λ⁴正是瑞利-金斯。2.2 MATLAB 里的稳定实现为什么用 expm1 而不是 exp(x)−1在 MATLAB 里实现普朗克函数最直觉的写法是1/(exp(x)-1)。但数值计算里当 x 很小时exp(x)-1会丢失有效数字。比如 x 1e-4exp(x) ≈ 1.000100005减去 1 之后只剩 5 位有效数字而 expm1(x) 是专门为这个场景设计的函数可以直接给出 1.000000005e-4。黑体辐射的长波端恰好就是 x 很小的区域。更省心的是当 x 很大时 expm1 溢出为 Inf倒数变成 0物理上对应短波端辐射强度趋近于零行为也正确。function M planck_lambda(lambda_um, T) % 普朗克光谱辐出度波长单位微米温度单位开尔文 % 输出单位: W/(m^2·um) h 6.626e-34; % 普朗克常数, J·s c 2.998e8; % 光速, m/s k 1.381e-23; % 玻尔兹曼常数, J/K lam lambda_um * 1e-6; % 微米换算成米 x h * c ./ (lam .* T); % 无量纲指数, 约等于 0.014388/(lam*T) M_per_m 2 * pi * h * c^2 ./ lam.^5 ./ expm1(x); M M_per_m * 1e-6; % W/m^3 换成 W/(m^2·um) end这段代码有两个值得注意的设计。第一lambda_um 传列向量、T 传行向量时lam .* T会借助隐式扩展生成一个 n×m 矩阵每一列对应一个温度后面画多曲线图时一次调用就能得到全部数据。第二2 * pi * h * c^2这个系数约等于 3.741e-16 W·m²是普朗克公式里最容易记错的部分合写避免手输常量出错。单位换算是另一个隐蔽的坑M_per_m 是按“每米波长”计算的 W/m³要显示成“每微米波长”必须乘 1e-6否则曲线整体偏离 6 个数量级。2.3 极限检验两条近似曲线如何夹住普朗克曲线验证一个物理公式的实现是否正确最直接的方法是看极限行为。把普朗克函数和瑞利-金斯、维恩近似画在同一张双对数图里T 1500; lam logspace(-1, 2, 500); % 0.1 到 100 微米 M_planck planck_lambda(lam, T); lam_m lam * 1e-6; % 转成米 M_rj (2 * pi * c * k * T ./ lam_m.^4) * 1e-6; % 瑞利-金斯 M_wien (2 * pi * h * c^2 ./ lam_m.^5 .* exp(-h*c./(lam_m*k*T))) * 1e-6; loglog(lam, M_planck, k-, LineWidth, 2); hold on; loglog(lam, M_rj, r--, LineWidth, 1.2); loglog(lam, M_wien, b-., LineWidth, 1.2); xlabel(波长 \lambda (μm)); ylabel(光谱辐出度 (W·m^{-2}·μm^{-1})); legend(普朗克, 瑞利-金斯, 维恩近似);双对数坐标下瑞利-金斯那条 λ⁻⁴ 曲线是一条斜率为 −4 的直线向左延伸到短波端时数值无限增大这就是紫外灾变的图像表达。维恩近似在长波端明显偏离、在短波端快速下坠。普朗克曲线在长波端贴着瑞利-金斯、在短波端贴着维恩中间平滑过渡。如果实现有问题最明显的症状就是长波端不再保持直线斜率 −4或者短波端出现异常的波动。3. 温度如何重塑曲线维恩位移与斯特藩-玻尔兹曼的数值验证3.1 峰值反比于温度维恩位移的推导与数值求解普朗克公式里温度只出现在指数项的分母上但观察曲线会发现温度一变峰的横坐标也在动。这个规律不是靠“看图”看出来的而是对 M_λ 求导得到的。令 x h c / (λ k T)峰值条件 d M_λ / dλ 0 可以整理成一个与温度无关的超越方程e^{-x} x/5 − 1 0这个方程里没有 T说明峰值对应的 x 是常数。用 fzero 解它x_peak fzero((x) exp(-x) x/5 - 1, 5); % 约 4.9651 lambda_max_um 1e6 * (h*c) / (k * x_peak * T);x_peak 约等于 4.9651于是 λ_max · T h c / (k x_peak) ≈ 2898 μm·K。这就是维恩位移定律温度翻倍峰值波长减半。T 300K 时峰值在 9.7μm 的远红外区室温物体辐射基本看不见T 5800K 时峰值在 0.5μm 附近正好落在可见光绿蓝区间——太阳光谱峰值在可见光附近不是巧合而是表面温度决定的。3.2 总量按 T⁴ 增长斯特藩-玻尔兹曼的数值验证把普朗克函数对整个波长积分得到总辐射出度 M_total σ T⁴σ ≈ 5.670e-8 W/(m²·K⁴)。这个四次方关系是“温度升高一点点辐射总量剧烈增长”的来源。用数值方法验证它比直接用公式更有手感——因为你会亲眼看到积分网格选择对结果的影响。lam_log logspace(-3, 6, 1e4); % 0.001 到 1e6 微米 M_at_300 planck_lambda(lam_log, 300); lam_m lam_log * 1e-6; sigma_num trapz(log(lam_m), M_at_300 .* lam_m); sigma_theory 5.670e-8;这里的关键是换元积分。普朗克函数输出的是“每微米”的辐出度积分时必须带着微米转米的因子 1e-6。直接trapz(lam_log, M_at_300)不只网格不理想单位也会混乱。换成对数网格后log(lam_m)是等距的自然对数轴被积函数变成 M_λ · λ梯形积分精度高而且长波拖尾被对数坐标压紧不会漏面积。运行后 sigma_num 和理论值对得上误差通常在 0.1% 以内。3.3 多曲线图把六个温度放进同一张图的正确姿势多曲线图是理解普朗克定律最直观的方式但线性坐标下它几乎没法看。先看实现再解释坐标选择T_list [300 500 1000 2000 3000 5800]; lam logspace(-1, 2.7, 2000); % 0.1 到 500 微米, 列向量 M planck_lambda(lam, T_list); % 2000x6 矩阵, 每列一条曲线 loglog(lam, M, LineWidth, 1.5); legend(cellstr(num2str(T_list)), Location, northeast); xlabel(波长 \lambda (μm)); ylabel(光谱辐出度 (W·m^{-2}·μm^{-1}));T_list 被传成行向量配合 lam 的列向量一次调用就把所有温度算完了。曲线在双对数坐标下呈现出清晰的“平移”趋势温度越高整条曲线越靠左上。这个视觉印象背后是普朗克公式的标度性质——M_λ / T⁵ 只是 λT 的函数所以把横轴换成 λT、纵轴换成 M_λ / T⁵六条曲线会完全重合。这个性质比公式本身更有价值它是后面一切不均匀刻度作图技巧的物理根据。4. 不均匀刻度图为什么线性坐标救不了多曲线图4.1 线性坐标下曲线“消失”的数学原因拿 300K 和 6000K 两条曲线对比。300K 的峰值约 3.1e-4 W/(m²·μm)6000K 的峰值约 9.9e7 W/(m²·μm)两者相差十一个数量级。线性纵轴把 6000K 的曲线画满整个坐标框时300K 的最高点还不到一个像素的高度看起来就是一条贴住零轴的水平线。这不是 MATLAB 绘图精度的问题而是物理量本身的动态范围太宽。横轴的情况类似。300K 的峰在 10μm 附近6000K 的峰在 0.48μm 附近而黑体辐射的有效范围从 0.1μm 一直延伸到 100μm 开外跨度超过三个数量级。线性横轴就像把北京到广州的距离摊在一把尺子上然后还要求你在上面标出一厘米的刻度——信息全被压扁了。解决思路不是放大局部而是换一把“不均匀的尺子”对数轴。4.2 半对数与双对数不均匀刻度的正确打开方式MATLAB 里切换坐标轴刻度是一行命令的事semilogx只把横轴变对数适合纵轴范围没那么夸张的场景loglog把两个轴都变对数适合黑体辐射这种横纵轴都跨数量级的场景。选择的标准很简单哪个物理量是幂律变化的哪个轴就用对数。双对数坐标还有个额外好处幂律关系会变成直线。瑞利-金斯区是 λ⁻⁴所以低温曲线的长波端在 loglog 图里是一条斜率为 −4 的直线维恩区近似指数衰减表现为快速下坠的弯曲段。你可以在图里检查任意一条曲线的最右端是否满足这个斜率如果偏差明显说明波长范围没给够或者数值实现有误。对数轴本质上就是一种不均匀刻度它把三个数量级的波长均匀地铺在屏幕宽度上这正是处理普朗克曲线的第一选择。4.3 别忘了乘 λ对数轴上的“面积守恒”很多人把横轴改成对数后直接拿 M_λ 当纵轴去分析“哪个波段贡献最多能量”这会得出错误的结论。原因在于对数横轴上曲线下的面积是 ∫ M_λ d(lnλ)而真实能量是 ∫ M_λ dλ。这两个积分相差一个 λ 的权重。正确做法是把纵轴改成 λ · M_λ这样 M_λ · λ · d(lnλ) M_λ · dλ面积才和能量对应semilogx(lam, M .* lam, LineWidth, 1.5);这个变换不只是绘图技巧它会改变峰的视觉位置。M_λ 的峰在 x 4.9651而 λ·M_λ 的峰在 x 3.9207——对 5800K 的黑体前者在 0.5μm 附近绿色后者在 0.63μm 附近橙红。肉眼觉得太阳“偏黄”和能量谱峰值偏长波是同一个现象的两面。做遥感或照明设计时如果分不清自己看的是哪种谱峰值位置会差出 20% 以上。4.4 可见光区间怎么突出不需要自定义分段坐标轴有人会想能不能做一个“0.4 到 0.7μm 线性放大、0.7 到 20μm 对数压缩”的分段坐标轴技术上可以思路是把波长映射到屏幕坐标对 λ 做分段插值线性段用线性映射远红段用 log10 映射然后用set(gca, XTick, ...)重新标刻度。但这套方案要手动解决网格线、刻度标签、轴范围三层问题属于典型的“为了一个视觉效果折腾一整天”。更务实的做法是对数坐标加可见光区间标注。0.38 到 0.78μm 恰好落在一个数量级内用对数轴看已经足够宽画上补丁区域就能明确位置yl ylim; patch([0.38 0.78 0.78 0.38], [yl(1) yl(1) yl(2) yl(2)], ... [0.9 0.9 0.8], FaceAlpha, 0.3);对数轴天然完成了“不均匀刻度”的职责短波被拉伸长波被压缩可见光区正好占据舒适的视觉宽度。除非期刊排版硬性规定坐标轴必须是线性分段样式否则对数轴加补丁是最省力又最不容易出错的方案。5. 避坑指南多曲线图里的系统性翻车现场5.1 现象低温曲线画出来是一条贴零的直线第一次画完多曲线图300K 那条几乎看不见很多人第一反应是“数据太小了”。原因上面已经算过300K 峰值只有 6000K 的约 3e-11 倍线性纵轴下物理上就不可能看到它。解决方法是改用 loglog或者对每条曲线做峰归一化M_norm M ./ max(M); % max 默认取每列最大值, 隐式扩展逐列归一 plot(lam, M_norm);归一化之后六条曲线峰值都拉到 1但峰的位置移动仍然保留适合用来演示维恩位移定律。注意归一化只对非零曲线有意义如果某列因为数值问题全部是 NaNmax 会返回 NaN归一化结果直接崩先排查数据。5.2 现象短波端曲线变成锯齿甚至断线短波端是指数变化最剧烈的区域如果波长用linspace(0, 100, 2000)生成第一点是 λ 0此时普朗克函数分子分母都是 0计算出 NaN曲线直接断开。即便避开零点线性采样在短波区也太稀疏指数下降段会被画成一折一折的折线。解决方法是恒用 logspace 生成波长并让下界远离零比如logspace(-1, 2.7, 2000)。对数采样让短波一侧的采样点密度远高于长波指数段的折线效应基本消失。如果仍然出现锯齿检查数据里有没有 Inf——温度高且波长极短时expm1 溢出为 Inf倒数变 0这是正确行为而不是 bug。5.3 现象峰值位置和 2898 μm·K 对不上维恩位移常数 2898 只适用于“按波长”形式的光谱辐射强度。很多人用频率形式的普朗克公式求峰值得到一个不同的常数然后怀疑公式抄错。两者都正确差别来自坐标变换的雅可比因子dν 和 dλ 之间相差一个 λ²/c 的比例峰的位置被这个因子推着走。公式形式峰值条件峰值常数典型应用场景M_λ按波长hc/(λkT) 4.9651λ_max · T ≈ 2898 μm·K光学测温、遥感、LED 光谱M_ν按频率hν/(kT) 2.8214ν_max /T ≈ 5.879×10¹⁰ Hz/K量子光学、射电天文画图前先确认自己画的是 M_λ 还是 M_ν再决定用哪个常数做自检。把频率谱强行折到波长轴后找峰值找出来的位置既不是 2898 也不是直觉上的“正确峰”这是初学者最容易掉的坑。5.4 现象数值积分验证斯特藩-玻尔兹曼总是差几个百分点用trapz(lam, M)直接对线性波长积分结果普遍偏小。主要误差来源不是截断而是采样线性网格下峰值附近的窄峰只落在一两个点上梯形法低估峰高长波尾部的贡献反而因为 λ⁻⁴ 衰减得很快截断到 1000μm 损失可以忽略。修正办法是换元到对数网格再积分也就是 3.2 节那段代码。积分上限推到 1e6μm 后300K 的尾部贡献已经可以忽略再继续推上限只会增加采样点数量不会明显改善精度。如果结果还是差把采样点数加到 1e5 再跑一次。5.5 现象默认配色让温度越高颜色越“冷”MATLAB 默认的 parula 色图不是温度顺序画六条温度曲线时300K 可能是蓝色5800K 可能是黄色视觉顺序和物理顺序完全错位。人眼读图时天然倾向“蓝 低温、红 高温”默认配色会造成系统性误读。手动指定色序set(gca, ColorOrder, turbo(numel(T_list)));新版 MATLAB 提供的 turbo 色图从深蓝渐变到深红和温度直觉一致。图例顺序也要和 ColorOrder 对应否则颜色和温度对不上血泪经验。这个细节在汇报场合比任何理论推导都容易被别人一眼看出问题。6. 横轴选波长还是选频率一个绕不开的设计取舍普朗克公式的两种形式——按波长和按频率——描述的是同一个物理对象但峰值位置、曲线形状、应用场景全都不一样。这不是物理矛盾而是“密度”这个概念依赖坐标系的必然结果辐射能量在波长轴上均匀分布和在频率轴上均匀分布是两回事雅可比因子 λ²/c 把它们区分开。工程上测温、遥感几乎都用波长形式因为传感器按波长采样量子光学和天文里频率形式更自然因为光子能量直接正比于频率。理解了这一点再看那些“为什么两个普朗克公式峰值对不上”的争论本质上只是坐标系选择问题。这个取舍的背后是一个更通用的习惯作图的坐标轴选择不是装饰而是物理信息的编码方式。坐标轴选错了结论就会错。我从那以后每次用 MATLAB 画普朗克曲线都强制走一遍自检先用 fzero 算出当前温度的理论峰值波长再看图里峰的位置对不对然后确认横轴跨度覆盖峰值两侧至少各一个数量级避免曲线被截断成“半边脸”。这两步加起来只要十秒钟却帮我避免了好几次温度标定上的乌龙。光谱数据是物理里面最容易“画出来像一回事、算起来另一回事”的东西你拿到的资料里如果有绘图脚本跑之前也先看一眼它的横轴是线性还是 logspace、纵轴有没有乘 λ确认完再信那条曲线。希望帮到你。本文还有配套的精品资源点击获取