如何提取scipy.integrate.quad返回结果中的有效积分值

问题背景 #

使用scipy.integrate.quad计算多上限数组的积分时，函数默认返回二元元组：第一个元素为有效积分计算结果，第二个元素为积分的绝对误差估计。当前代码直接存储了完整元组，需要提取所有元组的第一位有效结果，剔除误差项。
原有实现代码如下：

from scipy.integrate import quad
from typing import List

def integrand(z,alpha,beta,gamma):
    return (2**(-0.5 + 1/(2.*(1/(1 + z))**(3*gamma)))*3**(-0.5 + 3/(2.*(1/(1 + z))**(3*gamma)))*23**(1/(1 + z))**(-3*gamma)*np.exp(1 - (1/(1 + z))**(-3*gamma) - alpha/(2.*beta) + alpha/(2.*(1/(1 + z))**(3*gamma)*beta)))**(-1)


def multi_integrate(alpha,beta,gamma,lb:float, ub_ls:List[float]) -> List:
    results = []
    # loop through all upper bounds
    for ub in ub_ls:
        results.append(quad(integrand,lb,ub,args=(alpha,beta,gamma)))
    return results

lb = 0
ub_ls =z

def Hz_th(alpha,beta,gamma,z):
    return multi_integrate(alpha,beta,gamma, lb,ub_ls)

运行后返回结果为元组列表，格式如下：

[(0.00014626812936348427, 1.6239024498819118e-18),
 (0.000150015526635659, 1.6655069172217707e-18),
 (0.00015073611990008352, 1.6735071095573097e-18),
 (0.00015635608256829605, 1.7359012290751435e-18),
 (0.0001746488221728651, 1.938991435999794e-18)...

实现方法 #

有两种可行的提取方式，优先选择第一种，在计算阶段直接过滤无效数据，减少内存占用：

• 方法1：在循环计算时直接解包quad返回的元组，用下划线承接不需要的误差项，仅存储有效积分值。
  注意原代码缺少numpy导入，且ub_ls = z的定义位置不合理，修正后完整代码如下：

  import numpy as np
  from scipy.integrate import quad
  from typing import List
  
  def integrand(z,alpha,beta,gamma):
      return (2**(-0.5 + 1/(2.*(1/(1 + z))**(3*gamma)))*3**(-0.5 + 3/(2.*(1/(1 + z))**(3*gamma)))*23**(1/(1 + z))**(-3*gamma)*np.exp(1 - (1/(1 + z))**(-3*gamma) - alpha/(2.*beta) + alpha/(2.*(1/(1 + z))**(3*gamma)*beta)))**(-1)
  
  
  def multi_integrate(alpha,beta,gamma,lb:float, ub_ls:List[float]) -> List:
      results = []
      for ub in ub_ls:
          # 解包返回元组，忽略误差项
          integral_res, _ = quad(integrand,lb,ub,args=(alpha,beta,gamma))
          results.append(integral_res)
      return results
  
  lb = 0
  
  def Hz_th(alpha,beta,gamma,z):
      return multi_integrate(alpha,beta,gamma, lb, z)
  

• 方法2：如果已经生成了元组格式的结果列表，可通过列表推导式批量提取第一位元素：

  # raw_results为已存储的元组列表
  valid_results = [res[0] for res in raw_results]
  

效果 #

处理后返回的结果为纯数值列表，格式如下：

[0.00014626812936348427,
 0.000150015526635659,
 0.00015073611990008352,
 0.00015635608256829605,
 0.0001746488221728651,
 ...]

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