python-3.x - Python不准确的曲线拟合

Question

Temp    k(T)
298   6.66E-63
300   1.48E-62
350   3.58E-55
400   1.25E-49
450   2.57E-45
500   7.30E-42
550   4.90E-39
600   1.12E-36
650   1.11E-34
700   5.72E-33
750   1.75E-31
800   3.49E-30
850   4.92E-29
900   5.17E-28
950   4.24E-25
1000  2.83E-26

以上是给定的动力学数据，我正在尝试拟合这些数据并绘制相同的图。

曲线拟合

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd
plt.style.use('ggplot')

#Generate data
df=pd.read_excel('py_curvefit.xlsx')
T=df.Temp  #xdata
def reacKine(T,A,n,Ea):
    return A*((T/298)**n)*np.exp(-Ea/(0.008314*T))
kt=df['k(T)']  #ydata
#rectifying an erroneous value      
kt[14]=4.24*10**(-27)  
popt,pcov=curve_fit(reacKine,T,kt)
A,n,Ea=popt
plt.plot(T,np.log(kt),'g-',label='given data')
plt.plot(T,np.log(reacKine(T,*popt)),'ro',label='fit')
plt.xlabel('Temperature [K]')
plt.ylabel('log of reaction coefficient')
plt.legend(loc='best')
plt.show()

它表示未找到函数的最佳参数。我该如何纠正这个问题。我希望看到一个合适的。是因为指数项吗？

Answer 1

这是一个敏感问题(涉及指数时很典型)。对于这样的问题，对参数有一个很好的初始猜测是很重要的。

如果您试验这些参数，您会发现 A 必须非常小。 curve_fit 对所有参数使用的默认初始猜测为 1，而 1 对于 A 来说太大了。如果我使用 1e-10 作为 A

的初始猜测

popt, pcov = curve_fit(reacKine, T, kt, p0=(1e-10, 1, 1))

我从 curve_fit 得到以下错误:

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 800.

所以让我们将 maxfev 增加到，比如说，2000:

popt, pcov = curve_fit(reacKine, T, kt, p0=(1e-10, 1, 1), maxfev=2000)

我遇到了同样的错误。当我将它增加到 100000 时，函数成功了。

这是一个脚本，其中包含对 curve_fit 的更新调用，后跟脚本生成的绘图。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


T = np.array([298, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800,
              850, 900, 950, 1000])
kt = np.array([6.66e-63, 1.48e-62, 3.58e-55, 1.25e-49, 2.57e-45, 7.30e-42,
               4.90e-39, 1.12e-36, 1.11e-34, 5.72e-33, 1.75e-31, 3.49e-30,
               4.92e-29, 5.17e-28, 4.24e-27, 2.83e-26])

def reacKine(T,A,n,Ea):
    return A*((T/298)**n)*np.exp(-Ea/(0.008314*T))

popt, pcov = curve_fit(reacKine, T, kt, p0=(1e-10, 1, 1), maxfev=100000)


plt.plot(T, kt, '.', label='data')
tt = np.linspace(T[0], T[-1], 160)
kk = reacKine(tt, *popt)
semilogy = True
if semilogy:
    plt.semilogy(tt, kk, 'k-', alpha=0.3, label='fit')
    results_xy = (700, 1e-45)
else:
    plt.plot(tt, kk, 'k-', alpha=0.3, label='fit')
    results_xy = (300, 1.5e-26)

plt.annotate(xy=results_xy,
            s=('Fit Results:\n  $A\,$  = %.4g\n  $n\,$  = %.4g\n  $E_{a}$ = %.4g' %
               tuple(popt)))
plt.xlabel('T')
plt.ylabel('k(T)')
plt.legend(framealpha=1, shadow=True)
plt.show()

---

附言@MNewville 可能会建议一种更好的方法来使用 lmfit .