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

Question

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)

Answer 1

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 函数的减法会导致精度显着降低。