Python计算黑体辐射时np.exp触发overflow溢出报错原因求解

相关代码与报错信息 #

问题代码 #

import numpy as np

def BB_RAD(wavelength, T, d_lambda):
    c = 3 * 1E8 * 1E9 
    h = 6.626 * 1E-34
    k_B = 1.380649 * 1E-23
    X = (h*c)/(wavelength * k_B * T)
    u_lambda = ((2 * (np.pi) * (c**2) * h)/(wavelength**5)) *(1/(np.exp(X) - 1))
    f = u_lambda * d_lambda
    return f

wavelength = np.linspace(1,10000, 1000)
d_lambda = wavelength[2] - wavelength[1]
T = 26 + 273
E_lambda = BB_RAD(wavelength, T, d_lambda)

告警信息 #

BB_Radiation.py:14: RuntimeWarning: overflow encountered in exp
   u_lambda = ((2 * (np.pi) * (c**2) * h)/(wavelength**5)) *(1/(np.exp(X)- 1))

告警产生原因 #

• 物理规律层面：设置的黑体温度为299K（室温），根据维恩位移定律，该温度下黑体辐射的峰值波长约为10μm（即10000nm），设置的波长范围最低到1nm，属于远短于峰值的极短波区间，该区间内室温黑体的辐射强度理论上趋近于0。
• 数值计算层面：代码中计算的无量纲参数X = (h*c)/(wavelength * k_B * T)，在波长极小时数值会超过709，而双精度浮点数可支持的exp函数最大输入值约为709，超过该阈值后np.exp(X)的计算结果会超出浮点数存储上限，触发溢出告警。实际上该场景下1/(np.exp(X)-1)的结果已经可以近似为0，告警本身不会影响长波段的有效计算结果。

可选修复方案 #

如果需要消除告警，可以对大X的场景做近似处理，修改辐射强度计算逻辑：

X = (h*c)/(wavelength * k_B * T)
# X大于700时直接将分母设为无穷大，对应辐射强度为0
exp_denominator = np.where(X < 700, np.exp(X) - 1, np.inf)
u_lambda = ((2 * np.pi * (c**2) * h)/(wavelength**5)) * (1 / exp_denominator)

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