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)

Yes, I did. The answers that I got matched my example above, so it does seem to work! Thanks — Scott B, Jul 04 '15 at 01:02

score 1 · Accepted Answer · answered May 19 '15 at 01:26

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)

Nico Schlömer · Answer 2 · 2017-02-21T11:39:43.340

Using linearly spaced points is not very efficient. A better idea would be to turn to quadrature schemes for quadrilaterals, e.g., here. For example:

import numpy
import quadpy

def f(x):
   return (x[0] - 1) / ((x[0] - 1)**2 + x[1]**2)

quad = numpy.array([
    [-1, -1],
    [+1, -1],
    [+1, +1],
    [-1, +1],
    ])

scheme = quadpy.quadrilateral.Stroud(5)

val = quadpy.quadrilateral.integrate(f, quad, scheme)

Python/Scipy - Integrate with Quad Along Axis

2 Answers2