使用 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

integrate.trapezoid

用于一维积分。对于 2D 积分，您需要分别对阵列的每个轴进行积分。

我通过对被积函数进行平方来更改您的函数，以获得非零结果。这将帮助我们确认代码是否按预期工作。由于我没有您的数据文件，因此我还需要生成数据。

import numpy as np
from scipy import integrate

nvz = 1024 # evaluate at more points to see that
nph = 1024 # result matches that of quadrature
g = 1

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

costh = np.linspace(-1, 1, nvz)[:, np.newaxis]
phi = np.linspace(-np.pi, np.pi, nph)
sinth = np.sqrt(1 - costh**2)
integrand = (sinth * np.cos(phi))**2 * g
int_phi = integrate.trapezoid(integrand, dx=dphi, axis=-1)
res = integrate.trapezoid(int_phi, dx=dtheta)  # 4.180608973917667

比较：

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

integrate.dblquad(integrand, -1, 1, -np.pi, np.pi)
# (4.1887902047863905, 2.305878948871502e-09)