Solving double integration in python using scipy.integrate
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?
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.