使用scipy.curve_fit拟合普朗克辐射律出现popt=1、pcov=NaN问题求助

核心问题：变量与参数对应错误 #

你的代码里把自变量（频率v）和拟合参数（温度T）完全搞反了。普朗克公式的物理意义是：亮度是频率v和温度T的函数，即B(v, T) = (2hv³/c²) / (exp(hv/(kT)) - 1)，正确的对应关系是：

• 自变量为频率v
• 需要拟合的参数是温度T

你的planck函数定义中，错误将x当作温度、v当作频率，但调用curve_fit时却把频率传给x，拟合参数变成了v，这完全违背物理逻辑，直接导致拟合发散、出现除零警告。

修正步骤 #

1. 重新定义拟合函数 #

修正参数与变量的对应关系，同时处理低频区域的数值计算，避免除零：

import numpy as np
from scipy.optimize import curve_fit

# 使用SI定义的精确物理常量
h = 6.62607015e-34  # 普朗克常数
c = 299792458       # 光速精确值
k = 1.380649e-23    # 玻尔兹曼常数

def planck(v, T):
    # 计算指数项 hv/(kT)
    exponent = h * v / (k * T)
    # 处理低频区域（exponent<<1），用泰勒近似避免除零
    mask = exponent < 1e-6
    result = np.empty_like(v)
    # 泰勒展开近似：1/(exp(x)-1) ≈ 1/x - 1/2，避免分母趋近于0
    result[mask] = (2 * h * v[mask]**3) / (c**2 * (exponent[mask] + exponent[mask]**2 / 2))
    # 高频区域正常计算
    result[~mask] = (2 * h * v[~mask]**3) / (c**2 * (np.exp(exponent[~mask]) - 1))
    return result

2. 正确调用拟合函数并设置初始参数 #

CMB（宇宙微波背景）的温度约为2.7K，给拟合参数设置合理初始值，帮助算法收敛：

# 提取数据集
v_data = cmb['frequency (Hz)'].values
B_data = cmb['brightness (J/s/m^2/sr/Hz)'].values

# 拟合，指定初始参数p0
popt, pcov = curve_fit(planck, v_data, B_data, p0=[2.7])

print(f"拟合温度T = {popt[0]:.4f} K")
print(f"参数协方差矩阵:\n{pcov}")

3. 解决除零警告的关键 #

当hv/(kT) << 1（低频区域），exp(hv/(kT)) ≈ 1 + hv/(kT)，直接计算会导致分母趋近于0。通过泰勒近似处理低频数据，既能避免除零错误，又能保证计算精度。

4. 可选：数值缩放优化 #

如果数据数值范围差异极大（比如频率跨度从1e9到1e12 Hz），可对自变量缩放（如转成GHz），减少数值溢出风险：

v_scaled = v_data / 1e9  # 转换为GHz单位

def planck_scaled(v_ghz, T):
    v = v_ghz * 1e9
    exponent = h * v / (k * T)
    mask = exponent < 1e-6
    result = np.empty_like(v_ghz)
    result[mask] = (2 * h * v[mask]**3) / (c**2 * (exponent[mask] + exponent[mask]**2 / 2))
    result[~mask] = (2 * h * v[mask]**3) / (c**2 * (np.exp(exponent[mask]) - 1))
    return result

popt_scaled, pcov_scaled = curve_fit(planck_scaled, v_scaled, B_data, p0=[2.7])

关于对数拟合失效的原因 #

之前对数拟合无效，本质也是变量与参数对应错误，加上原始计算中存在数值不稳定（如出现零或负数导致取对数报错），修正变量对应关系后再尝试对数拟合即可生效。

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