python - 极限数组上的scipy.integrate.quad

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

我遇到一种情况，我需要评估成千上万个不同a，b的函数积分并求和。这需要很长时间才能遍历。我试图只给a和b提供四元组数组，希望四元组会返回相应的数组，但这没有用。

这是一个代码，它说明了我正在尝试做的事情，既可以使用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)

最佳答案

exp(-x*x)的积分是error function的缩放版本，因此您可以使用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函数的减法会导致精度显着下降。

关于python - 极限数组上的scipy.integrate.quad，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/29035107/