Python实现二维数组的二重积分

---

一、任务介绍

   在做海浪谱反演时，需要对某个函数进行定积分操作。虽然该函数的表达式已得知，但利用python的quad（）函数根据解释式求积分时，编程没有实现。因此选择利用梯形求解法对具体数值进行积分。
   给定一个函数z=f(x,y),该函数解析式不得知，但以确定二维NxM的数组Z.shape=[N,N]的形式给出。

	Z=[[1,2,3,...,7],...,[2,5,6,...,8]]

x和y的数值也给出。x.shape=N，y.shape=M。

1.任务目标：

需要求得给定数组z对x和y的二重积分。

---

二、梯形数值分析

   在高等数学中，最原始的对函数进行定积分求解方法，就是将函数分解为无穷小的矩形，求每个小矩形的面积后，将每个小矩形面积进行相加。加和结果便为积分结果。
   在获得数组的具体值后，实际就是获得了函数的每个值。因此对一维函数f(x)来说，就是将每个点对应的数值f(x_n)乘以区间dx。对二维数组来说，就是先对所有x进行相加，然后再对所有y进行相加。
   代码如下。注意，sum（）函数是默认对第一个维度求和。x默认为第一维度，因此直接对第一维度的dx进行求和。

x=np.linspace(0.001,0.2,100)
y=np.linspace(-np.pi,np.pi,50)
f_c=f_fun(x,y)
dx=(0.2-0.001)/100
dy=(2*np.pi/50)
inte_x=sum(dx*f_c)
print('int_x',inte_x.shape)
int_y=sum(dy*inte_x)
print('int',inte_y)

最后结果与直接使用trapz()函数结果几乎一致。

三、integrate.trapz（）函数

直接使用integrate.trapz()函数。代码如下
先对x求积分和先对y求积分，最终结果是一致的。

from scipy import integrate
int=integrate.trapz(f_v_c,x)
print(int.shape)
int=integrate.trapz(int,y)
print('int',int)

总结

在知道函数的每个具体值，但不知道函数的具体解析式时，可以采用梯形数值分析的方法来进行求解积分。dx和dy的数值越小，即分的程度越细致，最终结果越准确。