scipy 对具有变量边界的数组进行积分

Question

我正在尝试在点列表上集成一个函数，并将整个数组传递给一个集成函数，以便对事物进行矢量化。对于初学者来说，调用 scipy.integrate.quad 太慢了，因为我有大约 10 000 000 个点需要积分。使用 scipy.integrate.romberg 可以更快地完成该技巧，几乎是瞬时的，而quad则很慢，因为您必须循环遍历它或对其进行矢量化。

我的函数相当复杂，但出于演示目的，假设我想将 x^2 从 a 积分到 b，但 x 是一个用于计算 x 的标量数组。例如

将 numpy 导入为 np

from scipy.integrate import quad, romberg

def integrand(x, y):   

    return x**2 + y**2


quad(integrand, 0, 10, args=(10) # this fails since y is not a scalar

romberg(integrand, 0, 10)  # y works here, giving the integral over
                           # the entire range 

但这仅适用于固定边界。有没有办法做类似的事情

z = np.arange(20,30)
romberg(integrand, 0, z)  # Fails since the function doesn't seem to
                          # support variable bounds

我看到的唯一方法是在 numpy 中重新实现算法本身并使用它，这样我就可以拥有变量边界。有什么函数支持这样的功能吗？还有 romb，您必须直接提供被积函数的值和 dx 区间，但这对于我的复杂函数来说太不精确（marcum Q 函数，找不到任何实现，这可能是另一种点它的方法).

Answer 1

尝试评估特殊函数时的最佳方法是编写一个函数，使用该函数的属性在所有参数范围内快速准确地评估它。单一方法不太可能为所有参数范围提供准确（甚至稳定）的结果。在许多情况下，直接评估积分（如本例所示）几乎肯定会失败。

话虽这么说，在多个范围内评估积分的一般问题可以通过将积分转化为微分方程并求解来解决。大致步骤是

1. 给定一个积分 I(t)，我假设它是函数 f(x) 从 0 到 t 的积分[这可以推广到任意下限]，将其写为微分方程 dI/dt = f( x).
2. 在 0 到 t 的某个时间范围内，针对某些初始条件（此处为 I(0)），使用

   scipy.integrate.odeint()

   求解该微分方程。该范围应包含所有感兴趣的限制。采样的精细程度取决于函数以及需要评估的准确程度。
3. 结果将是我们输入的 t 集合从 0 到 t 的积分值。我们可以使用插值将其变成“连续”函数。例如，使用样条线我们可以定义

   i = scipy.interpolate.InterpolatedUnivariateSpline(t,I)

   。
4. 分别给定数组

   b

   和

   a

   中的一组上限和下限，那么我们可以一次性将它们全部评估为

   res=i(b)-i(a)

   。

这种方法是否适用于您的情况需要您仔细研究您的参数范围。另请注意，Marcum Q 函数涉及半无限积分。原则上这不是问题，只需将积分变换为有限范围内的一即可。例如，考虑变换 x->1/x。无法保证此方法对于您的问题在数值上是稳定的。