lmfit配合integrate.quad拟合NMR超导Hebel-Slichter峰报错求助

错误原因 #

你遇到的报错和integrate.quad传参无关，根源是你对T1Textended函数使用了np.vectorize包装：numpy.vectorize返回的函数对象签名包含可变参数*args，而lmfit的Model初始化时需要解析函数的显式参数列表，识别到可变参数就会抛出不支持的错误。

解决方案 #

第一步：移除np.vectorize，手动改造函数支持数组输入 #

np.vectorize本身只是循环的语法糖，没有性能提升，直接在函数内部对输入的Temperature数组做循环处理即可，保证函数参数列表是显式的，lmfit可以正常识别。
修改后的函数示例：

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from lmfit import Model, Parameters

# 被积函数保持不变
def T1Tfunc(En, Temperature , Gamma0 , Nfactor, Gap2 , Tc):
    kB = 8.617E-5
    Delta0 = kB * Tc * Gap2 / 2
    Delta1 = Delta0 * np.tanh(((Tc / Temperature)-1) ** 0.5)
    Gamma1 = Gamma0 * (Temperature / Tc) ** Nfactor
    Enp = En + 8.974E-6
    EnB = En + Gamma1 * 1j
    EnBp = Enp + Gamma1 * 1j
    Ns = (EnB / np.sqrt(EnB * EnB - Delta1 * Delta1))
    Nsp = (EnBp / np.sqrt(EnBp * EnBp - Delta1 * Delta1))
    Ms = (Delta1 / np.sqrt(EnB * EnB - Delta1 * Delta1))
    Msp = (Delta1 / np.sqrt(EnBp * EnBp - Delta1 * Delta1))
    FE = 1/(1 + np.exp(En/(kB*Temperature)))
    FEp = 1/(1 + np.exp(Enp/(kB*Temperature)))
    func = (np.real(Ns)*np.real(Nsp)+np.real(Ms)*np.real(Msp))*FE*(1-FEp)
    return func

# 改造后的主函数，移除vectorize，手动支持数组输入
def T1Textended(Temperature , Gamma0 , Nfactor, Gap2 , Tc , Koringaa , Koringab):
    kB = 8.61728E-5
    res = []
    # 统一处理单值或数组输入
    for T in np.atleast_1d(Temperature):
        if T < (0.1 * Tc):
            T1Te = 0
        elif T < Tc:
            # 直接做积分，无需单独定义T1T函数
            I = integrate.quad(T1Tfunc, 0, 0.1, args=(T , Gamma0 , Nfactor , Gap2 , Tc))[0]
            T1Te = I*(2/(kB*T)) * (Koringaa + Koringab * T)
        else:
            T1Te = Koringaa + Koringab * T
        res.append(T1Te)
    return np.array(res)

第二步：修正参数创建方法 #

你原来的params = Model.Parameters()是错误用法，参数对象可以直接从Model实例生成，或者导入Parameters类创建：

# 导入数据
filename = 'Rb2CsC60.txt'
data = np.loadtxt(filename, delimiter=',')
datax = data[:, 0]
datay = data[:, 1]
dataerr = data[:, 2]

# 初始化Model，此时不会再报错
HSmodel = Model(T1Textended)
print(HSmodel.param_names, HSmodel.independent_vars) # 可以正常输出参数名和自变量

# 创建参数
params = HSmodel.make_params()
params.add('Tc', value=32.2, vary=False)
params.add('Gamma0', value=1E-3, vary=True)
params.add('Nfactor', value=1, vary=False)
params.add('Gap2', value=4.25, vary=True)
params.add('Koringaa', value=1, vary=True)
params.add('Koringab', value=0, vary=True)

# 执行拟合，带实验误差权重
result = HSmodel.fit(datay, params, Temperature=datax, weights=1/dataerr**2)

# 输出拟合结果
print(result.fit_report())

# 绘制拟合结果
result.plot()
plt.show()

额外优化建议 #

如果拟合速度较慢，可以适当缩小integrate.quad的积分上限，或者设置epsabs、epsrel参数降低积分精度，提升拟合效率。

内容的提问来源于stack exchange，提问作者Ross Colman