Scott B - 1 year ago 112

Python 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 Source

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)
```