Solving double integration in python using scipy.integrate

Question

I want to compute this integral:

I have a data file providing values of cos(theta), phi and g.

I am trying to solve it using the trapezoid method of scipy.integrate. But I am unsure if this is the correct way since it is a double integration and g depends on both cos_theta and phi.

The code is as follows :

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)

Can someone please suggest me a way to solve it correctly?

Answer 1

scipy.integrate.trapezoid is for 1D integrals. If you have a 2D integral, you'd need to integrate over each axis of your array separately.

Since I don't have your data file, I will also need to generate the data. In particular, I'll assume g = (costh + phi)**2, but it can be any function.

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

Compare against:

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)

Note that you could also use integrate.simpson as a drop-in replacement for integrate.trapezoid.