使用SciPy对二元函数以数组值为积分上限积分的问题咨询

现有代码的核心问题 #

• 入参定义错误：Xntot_integrated 不需要接收x参数，x是积分过程的内部哑变量，该函数仅需要传入积分上限z一个入参
• 被积函数逻辑错误：你已经将Xneq_x定义为(1+np.exp(x))**-1的计算结果，后续调用Xneq_x(x)属于将数值作为函数调用，会触发类型错误
• quad参数传递错误：scipy.integrate.quad要求被积函数的第一个入参为积分变量，其余固定参数需要通过args字段传入
• 调用规范错误：你导入的是scipy.integrate模块，调用quad时需要写全integrate.quad，或者单独导入quad

修正后可运行代码 #

import numpy as np
from scipy import integrate 

# 这里假设Q是已知常量，你代码里未给出定义，运行前请替换为实际取值
Q = 1  

def K(x):
    b = 0.252
    return b*(((4/(x**3))+(3/(x**2))+1/x) + (4/(x**3) + 1/(x**2))*np.exp(-x))

def Xntot_integrand(x, z):
    Xneq_x = (1+np.exp(x))**(-1)
    return Xneq_x**2 * np.exp(K(z) - K(x)) * np.exp(x)

def Xntot_integrated(z):
    # 取quad返回的第一个值为积分结果，第二个为误差估计值可忽略
    return integrate.quad(Xntot_integrand, 0, z, args=(z,))[0]

# 对积分函数向量化，支持直接输入数组形式的z
Xntot_integrated_vec = np.vectorize(Xntot_integrated)

T_narrow = np.linspace(1,0.01,100)
z_narrow = Q/T_narrow

final_integrated_Xneq = Xntot_integrated_vec(z_narrow)

方案选择建议 #

如果你的z数组规模在千级以内，直接用上述向量化积分的方案即可，计算精度更高；如果需要处理十万级以上的z值，可以先预计算几十到上百组z对应的积分结果，再用scipy.interpolate做插值计算，能大幅提升运行速度，两种方案都可以满足需求，按你的实际计算规模选择即可。

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