Scott B - 1 year ago 123
Python Question

# Python/Scipy - Integrate with Quad Along Axis

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)
print(res2)
``````

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)
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download