python求极限函数_用Python（Debye模型）拟合带参数极限的积分函数
我试图将电阻率与温度数据拟合到金属电阻率的Bloch-Gruneisen公式中：
正如你所看到的，有一个带参数极限的积分函数。我不知道如何实现一个算法来运行最小二乘拟合。
我想到了：import matplotlib.pyplot as plt
import numpy as np
import pylab as pl
import scipy as sp
from scipy.optimize import leastsq
#retrieve data from file
data = pl.loadtxt('salita.txt')
Temp = data[:, 0]
Res = data[:, 2]
def debye_func(p, T, r):
rho0, AD, TD = p
coeff = AD*np.power(T, 5)/np.power(TD, 4)
f = np.power(x^5)/np.power(np.sinh(x), 2) #function to integrate
err_debye = r - rho0 - coeff * #integral???
return err_debye
p0 = sp.array([0.0001 , 0.00001, 50])
plsq = leastsq(debye_func, p0, args=(Temp, Res))
print plsq
我该怎么写呢？在
编辑：我的代码：
^{pr2}$
现在我得到一个值错误：Traceback (most recent call last):
File "debye.py", line 24, in
plsq = leastsq(debye_func, p0, args=(Temp, Res))
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py", line 348, in leastsq
m = _check_func('leastsq', 'func', func, x0, args, n)[0]
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py", line 14, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
File "debye.py", line 19, in debye_func
err_debye = r - rho0 - coeff * quad(debye_integrand, 0, TD/(2*T))
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 247, in quad
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 296, in _quad
if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
我认为这意味着我至少提供了一个它不能接受的论点，但我不知道如何修改我的代码。在
编辑2：我用Maxima解析解出了我的方程，我得到import matplotlib.pyplot as plt
import numpy as np
import pylab as pl
import scipy as sp
from scipy.optimize import leastsq
from scipy.integrate import quad
from scipy.special import zetac
from mpmath import polylog
#retrieve data from file
data = pl.loadtxt('salita.txt')
Temp = data[:, 0]
Res = data[:, 2]
def debye_integrand(x):
return np.power(x, 5)/np.power(np.sinh(x), 2)
def debye_func(p, T, r, integral):
rho0, AD, TD = p
coeff = AD*np.power(T, 5)/np.power(TD, 4)
den = np.exp(TD/T) -1
m1 = 5*((TD/(2*T))**4)*np.log(np.exp(TD/(2*T)+1)*(np.exp(TD/T)-1)+120*polylog(5, np.exp(TD/(T))*(1-np.exp(TD/(2*T)))
m2 = 120*(TD/(2*T))*polylog(4, np.exp(TD/(2*T)))*(np.exp(np.exp(TD/T))-1)+60*((TD/(2*T))**2)*polylog(3, np.exp(TD/(2*T))*(1-np.exp((TD/(2*T)))
m3 = 20*((TD/(2*T))**3)*polylog(2, np.exp(TD/(2*T))*(np.exp(TD/T)-1)+120**polylog(5, -np.exp(TD/(2*T)))*(1-np.exp(TD/T))
m4 = 120*(TD/(2*T))*polylog(4, -np.exp(TD/(2*T)))*(np.exp(TD/T)-1)+60*((TD/(2*T))**2)*polylog(3, -np.exp(TD/(2*T)))*(1-np.exp(TD/T))
m5 = 20*((TD/(2*T))**3)*polylog(2, -np.exp(TD/(2*T)))*(np.exp(TD/T)-1) -2*((TD/(2*T))**5)*np.exp(TD/T)
m6 = 5*((TD/(2*T))**4)*np.log(1-np.exp(TD/(2*T))*(np.exp(TD/T)-1)
zeta = 15.0*zetac(5)/2
integral = (m1+m2+m3+m4+m5+m6)/den +zeta
err_debye = r - rho0 - coeff * integral
return err_debye
#initizalied with Einstein model fit
p0 = sp.array([0.00001 , 0.0000001, 70.0])
plsq = leastsq(debye_func, p0, args=(Temp, Res))
print plsq
它用m2表示SyntaxError: invalid syntax。我试着用数值的方法来做循环，但是没有成功。在
如果你想试试，我的.txt文件是here。第一列是温度，第三列是电阻率。在
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！