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

Question

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

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

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

现在，对于 i 的每个值，我需要对 x 中的许多点求积，以找到 y[i,arr ] over arr，其中 arr 是 x 的子域。

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

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

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

我建议的解决方案是:

1. 为 y[i,arr2] 创建一个新的插值对象，例如 y2，其中 arr2 是一个比 足够大的子域code>arr.
2. 创建一个新的linspace x2对应于arr和arr2的交集，然后使用现有的功能集成方法，如 scipy.integrate.quadrature 在 x2 上集成 y2。
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 个小时，如果您只需要执行一次，我认为这是可以的。但时间在很大程度上取决于您的数据。