与函数和采样数据的双重积分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) 

输出数组似乎朝着对角线造成电网被定向，同时更应返回一个数组，看起来像模糊输入数组。

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