各位大佬好!我是搞凝聚态物理的。我要利用最小二乘法拟合实验数据,就是这个对应的公式有点怪 \frac{dI}{dV} =-M \cdot sign(V) \cdot\int_{-\infty}^{\infty}N(E)\frac{df(E+eV)}{dV}dE , 其中: N(E)=Re(\frac{E-i\Gamma}{\sqrt{(E-i\Gamma)^2+\Delta^2}}) ,f=\frac{1}{exp(\frac{E}{kT})+1}. \frac{dI}{dV} 和V的关系是实验测得的。为了把M, \Gamma ,\Delta,T 四个参量拟合出来,我编写了如下代码:
import numpy as np
from scipy import integrate as integ
from scipy import constants as C
from scipy.optimize import leastsq
import pandas as pd
#准备工作
#定义常数
e = C.elementary_charge
k = C.k
#定义函数
def STS(V,p):
M,T,Gamma,Delta=p #待拟合的参数
#定义N(E)
def N(E):
return np.real((E - 1j * Gamma) / np.sqrt((E - 1j * Gamma) ** 2 - Delta ** 2))
# 定义函数dfdV(E, V)
def dfdv(E, V):
return -e*np.exp((E + V*e)/(T*k))/(T*k*(np.exp((E + V*e)/(T*k)) + 1)**2)
#准备拟合,定义Dynes方未积分的形式
def dSTSdE(E,V):
return N(E)*dfdv(E,V)
result, error = integ.quad(lambda E: dSTSdE(E,V), -np.inf, np.inf)
result = -np.sign(V)*result*M #采用符号函数
return result
def residuals(p, V,data):
"""
实验数据和拟合函数之间的差,p为拟合需要找到的系数
"""
return data - STS(V, p)
#导入数据
frame = pd.read_excel('数据.xlsx')
AA=np.array(frame)[:,0]
CC=np.array(frame)[:,2]
# 调用leastsq进行数据拟合
# p0为拟合参数的初始值# args为需要拟合的实验数据
p0=[3.1284, 0.49, 0.04*10**(-3)*e, 1.34*10**(-3)*e] # 第一次猜测的函数拟合参数
plsq = leastsq(residuals, p0, args=(CC, AA))
print(plsq[0]) #拟合结果 [T,Gamma,Delta]
import pylab as pl
pl.plot(AA, CC,"o", label=u"原始实验数据")
pl.plot(AA, STS(AA, plsq[0]), label=u"拟合数据")
pl.legend()
pl.show()
然后就显示了如下问题:
runfile('E:/STUDY/数据拟合/data processing.py', wdir='E:/STUDY/数据拟合')
e:\study\数据拟合\data processing.py:22: RuntimeWarning: overflow encountered in exp
return -e*np.exp((E + V*e)/(T*k))/(T*k*(np.exp((E + V*e)/(T*k)) + 1)**2)
e:\study\数据拟合\data processing.py:22: RuntimeWarning: invalid value encountered in divide
return -e*np.exp((E + V*e)/(T*k))/(T*k*(np.exp((E + V*e)/(T*k)) + 1)**2)
Traceback (most recent call last):
File C:\Program Files\Spyder\pkgs\spyder_kernels\py3compat.py:356 in compat_exec
exec(code, globals, locals)
File e:\study\数据拟合\data processing.py:45
plsq = leastsq(residuals, p0, args=(CC, AA))
File C:\Program Files\Spyder\pkgs\scipy\optimize\_minpack_py.py:413 in leastsq
shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
File C:\Program Files\Spyder\pkgs\scipy\optimize\_minpack_py.py:26 in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
File e:\study\数据拟合\data processing.py:35 in residuals
return data - STS(V, p)
File e:\study\数据拟合\data processing.py:26 in STS
result, error = integ.quad(lambda E: dSTSdE(E,V), -np.inf, np.inf)
File C:\Program Files\Spyder\pkgs\scipy\integrate\_quadpack_py.py:463 in quad
retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
File C:\Program Files\Spyder\pkgs\scipy\integrate\_quadpack_py.py:577 in _quad
return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
TypeError: only size-1 arrays can be converted to Python scalars
是根据python数据处理相关教科书编写的代码,不知道问题出在哪里,万望各位大佬出手相助!/(ㄒoㄒ)/~~
相似问题