为什么我的 Python 代码使用 scipy.curve_fit() 计算普朗克辐射定律会产生“popt=1”和“pcov=inf”错误？

Question

SciPy 的 Curvefit 中未估计协方差

这是我的数据集：

frequency (Hz)  brightness (ergs/s/cm^2/sr/Hz)  brightness (J/s/m^2/sr/Hz)
float64 float64 float64
34473577711.372055  7.029471536390586e-16   7.029471536390586e-19
42896956937.69582   1.0253178228238486e-15  1.0253178228238486e-18
51322332225.44733   1.3544045476166584e-15  1.3544045476166584e-18
60344529880.18272   1.6902073280174815e-15  1.6902073280174815e-18
68767909106.5062    2.0125779972022745e-15  2.0125779972022743e-18
77780126454.10146   2.3148004995630144e-15  2.3148004995630145e-18
... ... ...
489996752265.52826  3.201319839821188e-16   3.201319839821188e-19
506039097962.6759   2.5968748350997043e-16  2.596874835099704e-19
523273092332.3638   2.0595903864583913e-16  2.0595903864583912e-19
539918248580.7806   1.7237876060575648e-16  1.7237876060575649e-19
557158231134.7507   1.3879848256567381e-16  1.3879848256567383e-19
573803387383.1646   1.0521820452559118e-16  1.0521820452559118e-19
591049358121.42 9.178609330955852e-17   9.178609330955852e-20

我尝试使用 CurveFit 来拟合普朗克辐射定律：

import numpy as np
from scipy.optimize import curve_fit

h=6.626*10e-34
c=3*10e8
k=1.38*10e-23
const1=2*h/(c**2)
const2=h/k

def planck(x,v):
    return const1*(v**3)*(1/((np.exp(const2*v/x))-1))

popt,pcov= curve_fit(planck, cmb['frequency (Hz)'],cmb['brightness (J/s/m^2/sr/Hz)'])
print(popt, pcov)

警告：

/tmp/ipykernel_2500/4072287013.py:11: RuntimeWarning: divide by zero encountered in divide
  return const1*(v**3)*(1/((np.exp((const2)*v/x))-1))

我得到

popt=1

和

pcov=nan

。现在函数中的指数项相差几个数量级。并且某些值不允许以数学方式近似该定律。我尝试使用法律的对数形式，但这也不起作用。我该如何克服这个问题？

Answer 1

这里有很多问题，包括变量被交换，不必要地重新定义物理常数，以及表达式在数值上高度不稳定。您需要使用

exp1m

来代替：

import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import h, c, k
from scipy.optimize import curve_fit

freq, brightness_erg, brightness_j = np.array((
    (34473577711.372055, 7.0294715363905860e-16, 7.0294715363905860e-19),
    (42896956937.695820, 1.0253178228238486e-15, 1.0253178228238486e-18),
    (51322332225.447330, 1.3544045476166584e-15, 1.3544045476166584e-18),
    (60344529880.182720, 1.6902073280174815e-15, 1.6902073280174815e-18),
    (68767909106.506200, 2.0125779972022745e-15, 2.0125779972022743e-18),
    (77780126454.101460, 2.3148004995630144e-15, 2.3148004995630145e-18),
    (489996752265.52826, 3.2013198398211880e-16, 3.2013198398211880e-19),
    (506039097962.67590, 2.5968748350997043e-16, 2.5968748350997040e-19),
    (523273092332.36380, 2.0595903864583913e-16, 2.0595903864583912e-19),
    (539918248580.78060, 1.7237876060575648e-16, 1.7237876060575649e-19),
    (557158231134.75070, 1.3879848256567381e-16, 1.3879848256567383e-19),
    (573803387383.16460, 1.0521820452559118e-16, 1.0521820452559118e-19),
    (591049358121.42000, 9.1786093309558520e-17, 9.1786093309558520e-20),
)).T


def planck(v: np.ndarray, T: float) -> np.ndarray:
    return 2*h/c/c * v**3 / np.expm1(h*v/k/T)

guess = 2.5,

(T,), _ = curve_fit(
    f=planck, xdata=freq, ydata=brightness_j, p0=guess, method='lm',
    # bounds=(0.1, np.inf),
)
print('T =', T)

fig, ax = plt.subplots()
v = np.linspace(freq.min(), freq.max(), 500)
ax.scatter(freq, brightness_j, label='data')
ax.plot(v, planck(v, *guess), label='guess')
ax.plot(v, planck(v, T), label='fit')
ax.legend()
plt.show()