使用 SciPy 在矩形网格上集成 2D 样本答案

【问题标题】：integrating 2D samples on a rectangular grid using SciPy使用 SciPy 在矩形网格上集成 2D 样本
【发布时间】：2020-07-20 16:16:11
【问题描述】：

SciPy 具有三种对样本进行一维积分的方法（trapz、simps 和 romb）和一种对函数进行二维积分的方法（dblquad），但它似乎没有进行二维积分的方法超过样本——甚至是矩形网格上的样本。

我看到的最接近的是 scipy.interpolate.RectBivariateSpline.integral - 您可以从矩形网格上的数据创建一个 RectBivariateSpline，然后将其整合。但是，这并不是非常快。

我想要比矩形方法更准确的方法（即只是总结所有内容）。比如说，我可以使用 2D Simpson 规则，方法是创建一个具有正确权重的数组，将其乘以我想要积分的数组，然后对结果求和。

但是，如果那里已经有更好的东西，我不想重新发明轮子。有吗？

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

使用一维规则两次。

>>> from scipy.integrate import simps
>>> import numpy as np
>>> x = np.linspace(0, 1, 20)
>>> y = np.linspace(0, 1, 30)
>>> z = np.cos(x[:,None])**4 + np.sin(y)**2
>>> simps(simps(z, y), x)
0.85134099743259539
>>> import sympy
>>> xx, yy = sympy.symbols('x y')
>>> sympy.integrate(sympy.cos(xx)**4 + sympy.sin(yy)**2, (xx, 0, 1), (yy, 0, 1)).evalf()
0.851349922021627

【讨论】：

好电话。对于我的问题，我一直在对相同维度的数组进行积分，因此为 simpson 规则创建一个权重数组（称为数组 simp）并执行 sum(simp*z) 对我来说更快——因为我只需要定义 simp一次。尽管如此，双 simps 方法还是很值得了解的。
有没有办法在 3D 中做到这一点？？
这是真的，因为要积分的函数确实是f(x,y)=g(x)+h(y)，结果是正确的，因为积分区间是[0,1] x[0,1].

【解决方案2】：

trapz 可以通过以下方式在 2D 中完成。用示意图绘制一个点网格，

整个网格的积分等于小区域dS的积分之和。梯形法则将小矩形 dS 上的积分近似为面积 dS 乘以 dS 角上的函数值的平均值，即网格点：

∫ f(x,y) dS = (f1 + f2 + f3 + f4)/4

其中 f1、f2、f3、f4 是矩形 dS 角上的数组值。

观察每个内部网格点输入整个积分的公式四次，因为这对于四个矩形很常见。不在拐角的边上的每个点，进入两次是两个矩形的共同点，每个拐角点只进入一次。因此，积分是通过以下函数在 numpy 中计算的：

def double_Integral(xmin, xmax, ymin, ymax, nx, ny, A):

    dS = ((xmax-xmin)/(nx-1)) * ((ymax-ymin)/(ny-1))

    A_Internal = A[1:-1, 1:-1]

    # sides: up, down, left, right
    (A_u, A_d, A_l, A_r) = (A[0, 1:-1], A[-1, 1:-1], A[1:-1, 0], A[1:-1, -1])

    # corners
    (A_ul, A_ur, A_dl, A_dr) = (A[0, 0], A[0, -1], A[-1, 0], A[-1, -1])

    return dS * (np.sum(A_Internal)\
                + 0.5 * (np.sum(A_u) + np.sum(A_d) + np.sum(A_l) + np.sum(A_r))\
                + 0.25 * (A_ul + A_ur + A_dl + A_dr))

在 David GG 给出的函数上进行测试：

x_min,x_max,n_points_x = (0,1,50)
y_min,y_max,n_points_y = (0,5,50)
x = np.linspace(x_min,x_max,n_points_x)
y = np.linspace(y_min,y_max,n_points_y)

def F(x,y):
    return x**4 * y

zz = F(x.reshape(-1,1),y.reshape(1,-1))

print(double_Integral(x_min, x_max, y_min, y_max, n_points_x, n_points_y, zz))

2.5017353157550444

其他方法（Simpson、Romberg 等）可以类似推导出来。

【讨论】：

【解决方案3】：

如果你在一个矩形上处理一个真正的二维积分，你会得到这样的东西

>>> import numpy as np
>>> from scipy.integrate import simps
>>> x_min,x_max,n_points_x = (0,1,50)
>>> y_min,y_max,n_points_y = (0,5,50)
>>> x = np.linspace(x_min,x_max,n_points_x)
>>> y = np.linspace(y_min,y_max,n_points_y)
>>> def F(x,y):
>>>     return x**4 * y
# We reshape to use broadcasting
>>> zz = F(x.reshape(-1,1),y.reshape(1,-1))
>>> zz.shape 
(50,50)
# We first integrate over x and then over y
>>> simps([simps(zz_x,x) for zz_x in zz],y) 
2.50005233

你可以和真实的结果比较

【讨论】：