使用 scipy integration.quad 重复评估积分

Posted

技术标签:

【中文标题】使用 scipy integration.quad 重复评估积分【英文标题】:Repeated evaluations of integral using scipy integrate.quad 【发布时间】:2021-09-09 15:48:00 【问题描述】:

我有一个函数,我需要估计它在许多不同间隔上的定积分,其中一些端点彼此非常接近,

\int_a^b_k f(x) dx

我可以为每个间隔重新调用integrate.quad,但这似乎效率很低,因为integrate.quad 将为每次调用多次评估该函数。有什么方法可以重用 quad 已经评估的点来更快地对quad 进行后续调用?

(如果不是 integrate.quad,可能是另一个 scipy 函数或使用不同的公共 python 库?)

或者也许我应该缓存我之前的整数值,然后调用quad 来计算新点与之前最近计算的点之间的差异?

还要注意,所选端点并非都是预先知道的,它们是自适应选择的。

【问题讨论】:

quad 可以提供很多额外信息(您可以尽我所能阅读),但我看不到返回(或提供)实际集成点和价值观。但是评估你的功能有多昂贵?为新案例缓存和修改值可能与“原始”计算一样昂贵。 也许您可以使用不同的集成方法。据说quad 的一个优点是它选择了一个最佳(或至少是好的)评估点数(我还没有研究过细节)。另一方面,对于trapezoid 之类的东西,您提供xy 值的数组。根据功能和所需精度,可能需要比quad 更多的点,但可以使用numpy 数组方法“一次性”评估它们。 谢谢。评估可能相当昂贵:O(1000) 三角函数评估,我需要它实时。我发现quad 可以返回这些子区间的子区间和积分近似值。我不确定使用什么方法来计算子区间上的值(我验证它不仅仅是梯形),但它们确实与完整积分的值相加。也许我可以做一次完整的积分,保留这些点及其值,然后对它们之间的点进行自己的插值。 序列bₖ,是单调的吗? 不一定,但假设单调 b_k 的部分解决方案仍然有用。 【参考方案1】:

对于具有不同端点的区间,您几乎没有希望使用通用自适应积分方法。这样做的原因是,一般来说,如果你稍微改变一下限制,一个整数值就会发生很大的变化。你必须以某种方式提供你的函数在边界处不是“疯狂”的信息,给出估计等,我不知道这样的方法。

可以做的是将积分转换为相同的区间,例如[0, 1],通过变量替换。您最终需要在 [0, 1] 上集成许多(参数化)函数。

\int_0^1 f_a(x)

然后,您可以尝试对函数调用进行矢量化处理,即一次针对一个点或一组点评估所有函数。

quadpy(我的一个项目)支持矢量化求积,可以显着加快集成速度。以x ** a 为例,用于许多不同的指数:

import numpy as np
import quadpy

a = np.arange(1, 7)

def f(x):
    return np.power.outer(x, a).T


val, err = quadpy.quad(f, 0.0, 1.0)

print(val)
print(err)
[0.5        0.33333333 0.25       0.2        0.16666667 0.14285714]
[1.87519478e-20 2.58599737e-20 4.45666434e-20 5.19063171e-20
 5.45158477e-20 4.09060649e-20]

请注意,在f 中,计算是矢量化的,而不是遍历a。这大大加快了计算速度。

【讨论】:

看起来这与我需要的相似,但是一次评估一个参数的许多值,而不是在许多不同的间隔端点进行评估。同样不幸的是,我的用例是自适应的,这意味着后续评估的间隔终点可能取决于先前计算点的积分,因此无论如何我都无法对所有终点的评估进行矢量化。 也许你应该在你的问题中更准确地说明你的问题是什么。 你的意思是我正在集成的功能到底是什么?我需要使用几个因应用程序而异的功能来做到这一点,所以我不能具体说。我希望很清楚,我的一般问题是能够有效地评估端点可能不同的单个函数的多个定积分。 不,我是指您对问题的描述。你写“我有一个函数,我需要在许多不同的点估计它的积分”。 “在不同点进行估计”是什么意思? 谢谢,我编辑了原始描述以澄清:我有一个函数,我需要在许多不同的区间估计它的定积分,其中一些端点彼此非常接近。跨度>

以上是关于使用 scipy integration.quad 重复评估积分的主要内容,如果未能解决你的问题,请参考以下文章

将 CX_Freeze 与 Scipy 一起使用:scipy.special._ufuncs.py

使用 scipy.optimize.minimize 提前停止损失函数

无法使用 scipy.stats

为啥不使用 Scipy 的 FFT 代码中的结果与 Scipy FFT 不相似?

我正在使用 Python3,但我无法安装 scipy

SciPy和App Engine