如何在Matlab中通过无穷积分数值计算斯特藩-玻尔兹曼常数？

你的代码核心问题 #

• 未定义温度变量 T：必须先指定一个具体温度值（如室温300K），否则代码会直接报错。
• 循环效率低下：MATLAB中百万级步长的for循环运行速度极慢，向量化运算或内置积分函数是更优选择。
• 积分区间合理性：你选择的10nm到10μm区间确实覆盖了常见温度下Planck辐射的主要能量贡献区域，但可以通过验证确认是否足够。

改进方案1：向量化运算替代循环 #

通过一次性计算所有波长对应的辐射能量，彻底消除循环开销，大幅提升运行速度：

h = 6.626e-34; 
c = 2.998e8;
k = 1.381e-23; 
T = 300;  % 设定温度，示例为300K

lambdaMin = 10e-9;
lambdaMax = 10e-6; 
nSteps = 1000000; 

% 生成等间隔波长数组
lambda = linspace(lambdaMin, lambdaMax, nSteps);
dLambda = lambda(2) - lambda(1);

% 向量化计算所有波长对应的辐射能量
E = (2 * h * c^2) ./ (lambda.^5 .* (exp(h * c ./ (lambda * k * T)) - 1));
I = sum(E * dLambda);

% 理论值对比
sigma_theory = 5.670374419e-8;
I_theory = sigma_theory * T^4;

fprintf('数值积分结果：%.10f W/m²\n理论总辐射量：%.10f W/m²\n相对误差：%.4f%%\n', I, I_theory, abs(I-I_theory)/I_theory*100);

改进方案2：使用MATLAB内置积分函数 integral() #

integral() 会自适应调整积分步长，处理数值稳定性问题（如波长趋近于0时的溢出风险），且无需手动设定步长，精度和效率都更优：

h = 6.626e-34; 
c = 2.998e8;
k = 1.381e-23; 
T = 300;

% 定义Planck辐射函数句柄
planck_fun = @(lambda) (2 * h * c^2) ./ (lambda.^5 .* (exp(h * c ./ (lambda * k * T)) - 1));

% 直接积分从0到无穷大（函数在两端趋近于0，可安全计算）
I = integral(planck_fun, 0, Inf);

% 计算斯特藩-玻尔兹曼常数
sigma_calculated = I / T^4;
sigma_theory = 5.670374419e-8;

fprintf('计算得到的斯特藩-玻尔兹曼常数：%.12f W/(m²K⁴)\n理论值：%.12f W/(m²K⁴)\n相对误差：%.4f%%\n', sigma_calculated, sigma_theory, abs(sigma_calculated-sigma_theory)/sigma_theory*100);

内容的提问来源于stack exchange，提问作者Anon