Python中带可变积分限的积分函数曲线拟合问题求助

需要对含积分的函数使用curve_fit求解系数c1至c7，被积函数依赖多个数组变量（defor、stress、init_def），积分限也随样本变化（init_def到defor），积分形式为time=f(defor,stress)，需对defor进行积分。原代码运行报错，请求解决。

原代码 #

import numpy as np
from scipy import integrate
from scipy.optimize import curve_fit
from scipy.integrate import quad

time =np.array([0, 18, 24, 42, 48, 66, 72, 0, 4, 22, 28, 46, 52, 70])
defor =np.array([0.11, 0.62, 0.73, 0.91, 1.0, 1.17, 1.22, 0.15, 0.26, 0.4, 0.43, 0.51, 0.51, 0.58])
init_def =np.array([0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])
stress =np.array([10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0,7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5])
temp=943
e=2.71828  
        
arg=zip(defor,stress,init_def)
    
# subint_funct- enter subintegral function with с1-с7 uncertain coefficients      
def subint_funct(arg, c1,c2, c3,c4,c5,c6,c7):
    return ((c1**(-1))*(temp**c2)*(stress**(-c3))*e**((c5-c6*stress)/temp)*
       ((defor+1)**(-c3))*(defor**c4)*e**(-((c6*stress+c7)*defor/temp))  ) 
    
# integration- integration of subint_funct with lower and upper integration bounds init_def and defor     
def integration(arg,c1,c2,c3,c4,c5,c6,c7):
    return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]

# curve fitting, dependent data-time
      
parameter = sp.curve_fit(integration,
                     arg,
                     time)
parameter=parameter[0]
print (parameter)

报错信息 #

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [11], in <cell line: 25>()
     22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
     23     return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]
---> 25 parameter = sp.curve_fit(integration,
     26                          arg,
     27                          time)
     28 parameter=parameter[0]
     29 print (parameter)

File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:789, in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    787 # Remove full_output from kwargs, otherwise we're passing it in twice.
    788 return_full = kwargs.pop('full_output', False)
---> 789 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    790 popt, pcov, infodict, errmsg, ier = res
    791 ysize = len(infodict['fvec'])

File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:410, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    408 if not isinstance(args, tuple):
    409     args = (args,)
---> 410 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    411 m = shape[0]
    413 if n > m:

File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:24, in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     22 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     23                 output_shape=None):
---> 24     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     25     if (output_shape is not None) and (shape(res) != output_shape):
     26         if (output_shape[0] != 1):

File ~\anaconda3\lib\site-packages\scipy\optimize\minpack.py:485, in _wrap_func.<locals>.func_wrapped(params)
    484 def func_wrapped(params):
---> 485     return func(xdata, *params) - ydata

Input In [11], in integration(arg, c1, c2, c3, c4, c5, c6, c7)
     22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
---> 23     return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]

Input In [11], in <listcomp>(.0)
     22 def integration(arg,c1,c2,c3,c4,c5,c6,c7):
---> 23     return [quad(subint_funct,init_def,defor,args= (c1,c2,c3,c4,c5,c6,c7) )[0] for defor,stress,init_def in arg ]

File ~\anaconda3\lib\site-packages\scipy\integrate\quadpack.py:351, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    348 flip, a, b = b < a, min(a, b), max(a, b)
    350 if weight is None:
---> 351     retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
    352                    points)
    353 else:
    354     if points is not None:

File ~\anaconda3\lib\site-packages\scipy\integrate\quadpack.py:463, in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    461 if points is None:
    462     if infbounds == 0:
---> 463         return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    464     else:
    465         return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

TypeError: only size-1 arrays can be converted to Python scalars

错误原因及修正方案 #

错误点分析 #

1. 被积函数参数错误：quad会将积分变量作为第一个参数传入被积函数，但原subint_funct用arg作为第一个参数，且直接引用全局数组stress、defor，导致quad传入的标量积分变量和全局数组冲突，触发类型错误。
2. 函数参数顺序不符合curve_fit要求：curve_fit要求拟合函数的格式为f(xdata, *params)，其中xdata是所有输入变量的集合，params是待拟合的系数。原integration函数参数顺序颠倒，且用zip生成的迭代器只能被遍历一次，无法满足curve_fit多次调用的需求。
3. 未正确导入模块：原代码使用sp.curve_fit但未定义sp，应直接使用导入的curve_fit。

修正后的代码 #

import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit

# 数据准备
time = np.array([0, 18, 24, 42, 48, 66, 72, 0, 4, 22, 28, 46, 52, 70])
defor = np.array([0.11, 0.62, 0.73, 0.91, 1.0, 1.17, 1.22, 0.15, 0.26, 0.4, 0.43, 0.51, 0.51, 0.58])
init_def = np.array([0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])
stress = np.array([10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5])
temp = 943
e = np.e  # 直接用numpy的自然常数更准确

# 被积函数：第一个参数是积分变量x，后续是拟合系数和样本参数
def subint_funct(x, c1, c2, c3, c4, c5, c6, c7, stress_val, temp_val):
    return (
        (c1 ** (-1)) * (temp_val ** c2) * (stress_val ** (-c3)) * np.exp((c5 - c6 * stress_val) / temp_val)
        * ((x + 1) ** (-c3)) * (x ** c4) * np.exp(-((c6 * stress_val + c7) * x / temp_val))
    )

# 拟合函数：接收xdata（包含defor, stress, init_def）和拟合系数，返回每个样本的积分结果
def integration(xdata, c1, c2, c3, c4, c5, c6, c7):
    results = []
    # 遍历每个样本的参数
    for d, s, i_d in zip(xdata[0], xdata[1], xdata[2]):
        # 对当前样本从init_def到defor积分
        integral, _ = quad(subint_funct, i_d, d, args=(c1, c2, c3, c4, c5, c6, c7, s, temp))
        results.append(integral)
    return np.array(results)

# 将输入数据整理成curve_fit要求的格式：(defor数组, stress数组, init_def数组)
xdata = (defor, stress, init_def)

# 初始参数猜测（需根据实际情况调整，避免拟合不收敛）
p0 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

# 执行曲线拟合
popt, pcov = curve_fit(integration, xdata, time, p0=p0)

print("拟合得到的系数c1-c7：")
print(popt)

关键修正说明 #

1. 被积函数重构：将积分变量x作为第一个参数，样本参数stress_val、temp_val通过args传入，避免引用全局变量，确保quad能正确传入标量积分变量。
2. 拟合函数调整：按照curve_fit要求的参数顺序编写，将输入数据整理为元组形式，遍历每个样本执行积分，返回数组形式的结果（与time数组形状一致）。
3. 初始参数设置：添加p0参数为拟合提供初始猜测，提高拟合收敛的概率（可根据实际物理意义调整初始值）。
4. 使用np.e替代手动定义的e：保证自然常数的准确性。

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