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。第一列是温度，第三列是电阻率。在