Python中的函数和采样数据的双重积分

Question

我希望找到一种使用numpy trapz或scipy堆栈中类似的函数对采样数据进行双重积分的方法。特别地，我想计算以下函数：

其中f(x',y')是采样数组，F(x, y)是相同维度的数组。

这是我的尝试，但结果不正确：

def integrate_2D(f, x, y):
    def integral(f, x, y, x0, y0):
        F_i = np.trapz(np.trapz(np.arcsinh(1/np.sqrt((x-x0+0.01)**2+(y-y0+0.01)**2)) * f, x), y)
        return F_i
    sigma = 1.0
    F = [[integral(f, x, y, x0, y0) for x0 in x] for y0 in y]
    return F

xlin = np.linspace(0, 10, 100)
ylin = np.linspace(0, 10, 100)
X,Y = np.meshgrid(xlin, ylin)

f = 1.0 * np.exp(-((X - 8.)**2 + (Y - 8)**2)/(8.0))
f += 0.5 * np.exp(-((X - 1)**2 + (Y - 9)**2)/(10.0))

F = integrate_2D(f, xlin, ylin)

输出数组似乎是朝着结果网格的对角线方向定位的，而实际上应该返回一个看起来像模糊输入数组的数组。

- Gregzor

Answer 1

我能理解您的意图，但是嵌套可能会隐藏逻辑。尝试使用类似以下代码：

def int_2D( x, y, xlo=0.0, xhi=10.0, ylo=0.0, yhi=10.0, Nx=100, Ny=100 ):

    # construct f(x,y) for given limits
    #-----------------------------------
    xlin = np.linspace(xlo, xhi, Nx)
    ylin = np.linspace(ylo, yhi, Ny)
    X,Y = np.meshgrid(xlin, ylin)

    f = 1.0 * np.exp(-((X - 8.)**2 + (Y - 8)**2)/(8.0))
    f += 0.5 * np.exp(-((X - 1)**2 + (Y - 9)**2)/(10.0))

    # construct 2-D integrand
    #-----------------------------------
    m = np.sqrt( (x - X)**2 + (y - Y)**2 )
    y = 1.0 / np.arcsinh( m ) * f

    # do a 1-D integral over every row
    #-----------------------------------
    I = np.zeros( Ny )
    for i in range(Ny):
        I[i] = np.trapz( y[i,:], xlin )

    # then an integral over the result
    #-----------------------------------    
    F = np.trapz( I, ylin )

    return F

- Gabriel

---

是的！就是这样。我希望我能够添加一些图片来展示它的外观！ - Gregzor