How to integrate an array of functions?

Question

I'm trying to integrate each element of an array to output an array of the same size. Below, X and Y are meshgrid arrays, and trying to integrate returns the error "only size-1 arrays can be converted to Python scalars."

def integral_x(p_i, p_j):
    def integrand(s):
        return ((2*(X - (p_j.xa - math.sin(p_j.beta) * s))) /
                ((X - (p_j.xa - math.sin(p_j.beta) * s))**2 +
                 (Y - (p_j.ya + math.cos(p_j.beta) * s))**2))
    return integrate.quad(integrand, 0.0, p_j.length)[0]

I've tried using numpy.vectorize but it seems that integrate.quad just doesn't work with an array input. Is there any way to integrate an array besides using loops?

Edit: Quick example

N = 3                               # Number of points in each direction
x_start, x_end = -2.5, 2.5            # x-direction boundaries
y_start, y_end = -3.0, 3.0            # y-direction boundaries
x = np.linspace(x_start, x_end, N)    # 1D-array for x
y = np.linspace(y_start, y_end, N)    # 1D-array for y
X, Y = np.meshgrid(x, y)

def integral_x():
    def integrand(s):
        return ((2*(X - s))) /(((X - s))**2 +(Y - s)**2)
    return (integrate.quad((integrand), 0.0, 10)[0])

Answer 1

Just use integrate.quad_vec rather than integrate.quad. For example,

from scipy.integrate import quad_vec

N = 3                               # Number of points in each direction
x_start, x_end = -2.5, 2.5            # x-direction boundaries
y_start, y_end = -3.0, 3.0            # y-direction boundaries
x = np.linspace(x_start, x_end, N)    # 1D-array for x
y = np.linspace(y_start, y_end, N)    # 1D-array for y
X, Y = np.meshgrid(x, y)

def integral_x():
    def integrand(s):
        return ((2*(X - s))) /(((X - s))**2 +(Y - s)**2)
    return (quad_vec((integrand), 0.0, 10)[0])

print(integral_x())

results in the printed output

[[-1.45895275 -1.04304777 -0.1282991 ]
 [-2.53152698        -inf  0.60485952]
 [-2.686371   -3.58488316 -3.98234315]]