Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

huangapple go评论62阅读模式
英文:

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

问题

协方差未在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=1pcov=nan。现在函数中的指数项相差几个数量级。并且一些值不允许数学上近似这个定律。我尝试使用该定律的对数形式,但那也不起作用。我该如何克服这个问题?

英文:

Covariance not estimated in SciPy's Curvefit

Here's my dataset:

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

I tried to use CurveFit to fit this to Planck's Radiation Law:

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)

Warning:

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

I get popt=1 and pcov=nan. Now the exponential term in the function differs by several orders of magnitude. And some of the values don't permit to approximate the law mathematically. I tried using the logarithmic form of the law but that doesn't work either. How can I overcome this problem?

答案1

得分: 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()

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

英文:

A lot of problems here, including that your variables were swapped, you're needlessly redefining physical constants, and your expression was highly numerically unstable. You need to use exp1m instead:

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()

Why does my Python code using scipy.curve_fit() for Planck's Radiation Law produce 'popt=1' and 'pcov=inf' errors?

huangapple
  • 本文由 发表于 2023年6月1日 21:28:27
  • 转载请务必保留本文链接:https://go.coder-hub.com/76382415.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定