使用scipy.integrate在python中解决双重积分

Question

我想计算这个积分：

我有一个数据文件，里面提供了cos(theta)、phi和g的值。

我正在尝试使用scipy.integrate中的梯形法来解决这个问题。但我不确定这样做是否正确，因为这是一个双重积分，而且g是依赖于cos_theta和phi的。

代码如下：

nvz = 256
nph = 256
dOmega = (2.0/nvz) * (2*np.pi / nph)
dphi = (2*np.pi / nph)
dtheta = (2.0/nvz)
cos_theta = file[:,0]
sin_theta = np.sqrt(1-cos_theta**2)
phi = file[:,1]
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
g = file[:,2]

integrate.trapezoid(sin_theta*cos_phi*g, dx = dOmega)

有没有人能给我建议一个正确的解决方法？

Answer 1

scipy.integrate.trapezoid 是用来计算一维积分的。如果你要计算二维积分，就需要分别对数组的每个轴进行积分。

因为我没有你的数据文件，所以我也需要生成一些数据。具体来说，我会假设 g = (costh + phi)**2，但其实可以用任何函数。

import numpy as np
from scipy import integrate

nvz = 2560  # I'll evaluate at more points just to see closer
nph = 2560  # agreement with quadrature

dphi = (2*np.pi / nph)
dtheta = (2.0/nvz)

# align `costh` along axis 0 and `phi` along axis 1
costh = np.linspace(-1, 1, nvz)[:, np.newaxis]
phi = np.linspace(-np.pi, np.pi, nph)[np.newaxis, :]

# g can be any array broadcastable to shape `(nvz, nph)`
g = (phi + costh)**2

# integrate along axis 1, then integrate along remaining axis
sinth = np.sqrt(1 - costh**2)
integrand = sinth * np.cos(phi) * g
int_phi = integrate.trapezoid(integrand, dx=dphi, axis=1)
res = integrate.trapezoid(int_phi, dx=dtheta, axis=0)
res  # -19.72363915277977

可以和以下内容进行比较：

def integrand(phi, costh):
    sinth = np.sqrt(1 - costh**2)
    g = (phi + costh)**2
    return sinth * np.cos(phi) * g

integrate.dblquad(integrand, -1, 1, -np.pi, np.pi)
# (-19.73920880217874, 1.2569373097903735e-08)

注意，你也可以用 integrate.simpson 来替代 integrate.trapezoid，效果是一样的。