python - 数据功能在其域的一小部分上的多重集成 - 准确性和效率

Question

作为我研究的一部分，我需要将数据定义的函数多次集成到其域的小子集上。这是一篇很长的文章:以下三个问题中任何一个的答案都将被确认!

例如，假设我有一个很大的一维域 x ，和一个函数 y集成到 x 的某个子集上。

在一个典型示例中，域 x将有 1000 个网格点，其中步骤 ' i '我的数据函数y[i,:]是一个 numpy 数组，大小与 x 相同。通常y将是形状为 (1000,1000) 的 numpy 数组.

现在，对于 i 的每个值，我需要对 x 中的许多点采用正交。 ，求 y[i,arr] 的积分超过arr ，其中arr是 x 的子域。

这是我的第一个问题:当arr时很小(比如说3点)，像scipy.integrate.cumtrapz这样的方法不会给出一个好的近似值 - y[i,arr] 中只有三个值。

每一步i ，必须对 arr 进行此类集成范围从大小 3 到大约大小 200。这些集成的输出用于更新 y[i+1,:] ，所以我相信由于我当前使用 cumtrapz 而引入了很多错误。 .

编辑:非常感谢@Fabian Rost，他为问题 1 提供了答案:是否确实引入了错误。他还建议使用线性插值(如下面的问题 2 所示)，并估计这种技术需要多长时间。我想真正剩下的问题是是否有比提议的技术更快的技术。

我建议的解决方案是:

1. 创建一个新的插值对象，例如 y2 ，对于y[i,arr2] ，其中arr2是一个比 arr 足够大的子域.
2. 创建一个新的linspace x2对应于 arr 的交集与 arr2 ，然后使用现有的函数集成方法，如 scipy.integrate.quadrature整合y2超过x2 .
3. 第 2 步的结果可能非常接近 y[i,arr] 的积分.

问题 2:所有这些步骤都是必要的吗？也就是说，是否有一个内置程序可以为我完成这一切？

问题 3:我相信如果我想避免错误，我必须进行这些插值->积分，超过 1000 次迭代，每次迭代至少 200 个子域。这显然会变得相当昂贵。有没有更 Pythonic 的方法来做到这一点？

非常感谢您回答这些问题!非常感谢您的阅读。

Answer 1

假设线性插值是一个很好的模型，您可以定义一个连续数据函数并使用 scipy.integrate.quad 进行积分，如下所示:

import scipy as sp
import scipy.integrate

xp = sp.linspace(0, 1000, 1000)
yp = sp.randn(1000)

datafunc = lambda x: sp.interp(x, xp, yp)

sp.integrate.quad(datafunc, 3, 1000)

根据域大小，集成在我的计算机上需要 2 到 4 毫秒。这意味着 1000 * 200 集成需要大约 4 小时，如果您只需要执行一次，我认为这是可以的。但时间在很大程度上取决于您的数据。