python - Quadpy中具有非常数限制的双积分？

Question

我需要集成一个在 2Dims （z和 radius r）中定义的函数，对此我没有表达式。

我可以在任何位置查询函数(z,r)并获取返回值。

我有跨 z: 的积分范围[-z_range, z_range]，我将其分为以下N_z几点：

z_is = -z_range + (np.arange(N_z) + 0.5) * (2.*z_range/N_z)

对于 中的每个值z_is，要整合的范围r是：[0, r_thresh_at_this_z]。

r_thresh_at_this_z是从被z调用的值中获得的z_i：

def get_r_thresh(z_i):
    return expression_of_z_i # returns a positive float.

所以径向积分的范围取决于 上的值z。

我的函数 f 为：

def f(r,z):
    return interpolator_for_f(r,z)

我想以quadpy最有效的方式使用这个包，因为它被创建为能够以这种方式使用。

我正在考虑使用 for 循环遍历，并对 中的每个值z_is执行高斯正交。[0, r_thres_for_that_z]z_is

我可以使用：

results = np.zeros((N_z))
errs = np.zeros((N_z))
for i in range(N_z):
   def f(r):
      return interpolator_for_f(r,z_is[i])
   results[i], errs[i] = quadpy.quad(f, [0, r_thresh_at_this_z])

但我觉得 for 循环并不是使用 quadpy 的最有效方式。

你能告诉我我在只用快速的 numpy 数组做这个积分时缺少什么，所以没有 for 循环？

[我已经阅读x了函数输入的形状f。在我的情况下d = 1，因为它是跨 r 的线积分。 n = N_z因为我想执行N_z这样的线积分，然后我将它们相加以获得 1 个单标量，即（整个）双积分的结果。

p = 1000因为假设我想要 1000 个积分点r，对于z.

所以我需要在N_z * 1000点处对函数进行采样。

函数f应返回一个数组形状(N_z, 1000)

这些参数的识别是否有任何帮助？]

谢谢！

标签： pythonintegrationnumerical-methodsnumerical-integration

Answer 1

Quadpy 作者在这里。

我想以最有效的方式使用 quadpy 包，因为它被创建为能够以这种方式使用。

事实并非如此。Quadpy 有助于在各种规范域中集成功能，但在任意二维域上集成的功能尚不存在。doublequad/nquad 应该在 quadpy 中，但还没有。