scipy.integrate.quad 关于极限数组

Question

quad 来自 scipy.integrate 需要参数 func、a、b。其中 func 是要积分的函数，a 和 b 分别是积分下限和积分上限。 a 和 b 必须是数字。

我有一种情况需要计算一个函数对数十万个不同的 a、b 的积分并对结果求和。这需要很长时间才能循环。我试图只为 a 和 b 提供四边形数组，希望四边形可以 return 相应的数组，但这没有用。

这是一个代码，说明了我正在尝试做的事情，Python 循环有效，但速度很慢，而我尝试进行矢量化却无效。关于如何快速（numpy-is）解决这个问题的任何建议？

import numpy as np
from scipy.integrate import quad

# The function I need to integrate:
def f(x):
    return np.exp(-x*x)

# The large lists of different limits:
a_list = np.linspace(1, 100, 1e5)
b_list = np.linspace(2, 200, 1e5)

# Slow loop:
total_slow = 0  # Sum of all the integrals.
for a, b in zip(a_list, b_list):
    total_slow += quad(f, a, b)[0]
        # (quad returns a tuple where the first index is the result,
        # therefore the [0])

# Vectorized approach (which doesn't work):
total_fast = np.sum(quad(f, a_list, b_list)[0])

"""This error is returned:

 line 329, in _quad
    if (b != Inf and a != -Inf):
ValueError: The truth value of an array with more than 
one element is ambiguous. Use a.any() or a.all()
"""

编辑：

我需要集成的实际功能（rho）还包含另外两个因素。 rho0是与a_list、b_list等长的数组。 H 是一个标量。

def rho(x, rho0, H):
    return rho0 * np.exp(- x*x / (2*H*H))

编辑2：

分析不同的解决方案。 ´space_sylinder` 是发生积分的函数。 Warren Weckesser 的建议与通过简单的分析函数传递数组一样快，比慢 Python 循环快约 500 倍（注意程序甚至没有完成的调用次数，它仍然使用了 657 秒).

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
# Analytic (but wrong) approximation to the solution of the integral:
      108    1.850    0.017    2.583    0.024 DensityMap.py:473(space_sylinder)
# Slow python loop using scipy.integrate.quad:
       69   19.223    0.279  657.647    9.531 DensityMap.py:474(space_sylinder)
# Vectorized scipy.special.erf (Warren Weckesser's suggestion):
      108    1.786    0.017    2.517    0.023 DensityMap.py:475(space_sylinder)

Answer 1

exp(-x*x) 的积分是 error function, so you can use scipy.special.erf 的缩放版本，用于计算积分。给定标量 a 和 b，函数从 a 到 b 的积分是 0.5*np.sqrt(np.pi)*(erf(b) - erf(a)).

erf 是一个 "ufunc"，这意味着它处理数组参数。给定 a_list 和 b_list，您的计算可以写成

total = 0.5*np.sqrt(np.pi)*(erf(b_list) - erf(a_list)).sum()

函数 rho 也可以用 erf 处理，方法是使用适当的缩放比例：

g = np.sqrt(2)*H
total = g*rho0*0.5*np.sqrt(np.pi)*(erf(b_list/g) - erf(a_list/g)).sum()

在依赖它之前，对照你的慢速解决方案检查它。对于某些值，erf 函数的减法会导致精度显着下降。

scipy.integrate.quad 关于极限数组

scipy.integrate.quad on arrays of limits

python

arrays

numpy

scipy

integral