Python/Scipy - Integrate with Quad Along Axis

Question

I have a 2D array. The "xy" plane is a grid from (-1,-1) to (1,1). I want to compute and integral at each point where the function depends on the coordinates of the point.

I know that with discrete data I can use simps or trapz and specify an axis to integrate along (see example). Is this possible with scipy.integrate.quad without using an ugly loop as shown below?

import numpy as np
import scipy as sp
import scipy.integrate

x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X,Y = np.meshgrid(x, y)

z = np.linspace(1, 10, 100)

# Integrate discrete values using simps
def func(z):
   return (X - z) / ((X - z)**2. + Y**2)
res1 = sp.integrate.simps(func(z.reshape(-1, 1, 1)), z, axis=0)
print(res1)

# Integrate the function using quad at each point in the xy plane
res2 = np.zeros(X.shape)
for i in range(res2.shape[0]):
   for j in range(res2.shape[1]):
      def func2(z):
         return (X[i,j] - z) / ((X[i,j] - z)**2. + Y[i,j]**2)
      res2[i,j] = sp.integrate.quad(func2, 1, 10)[0]
print(res2)

Answer 1

Using the Cubature method from Prof. Steven Johnson, which I wrapped using Cython, you can achieve the integration at once doing:

import numpy as np

from cubature import cubature

x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X,Y = np.meshgrid(x, y)

z = np.linspace(1, 10, 100)

def func(z):
   return (X.ravel() - z) / ((X.ravel() - z)**2. + Y.ravel()**2)

res = cubature(1, func, np.array([1.]), np.array([10.]))[0].reshape(X.shape)