ENH: scipy.integrate: vectorize quad

[cdeil]

The scipy.interpolate functions can be evaluated on grids of points.
In contrast, the scipy.integrate functions currently can only be evaluated for a single bin.

Would it be possible to change the scipy.integrate functions to accepts grids of bins as input and return grids of bins as output?

This is commonly needed when discretising models onto a grid (see e.g. here or here).

Assuming this is desired functionality for scipy.integrate, is it possible to extend the current functions in a backwards-compatible way?

[mdhaber]

For quad, it seems reasonable to accept array_like limits of integration a and b for user convenience. I think we need to do this in stats.rv_continuous when only _pdf is defined. Initially, I would not attempt to be fancy and break up the integral into pieces and combine them; I would simply loop or use np.vectorize. With appropriate a and b, this could be used either to calculate the integral within bins or do a cumulative integral. Should be no problem with backward compatibility, since limits of integration can't be arrays right now, and only array input would change the output.

[mdhaber]

@Kai-Striega Here's another one that seems really simple - just vectorize for array a and b (which can be broadcasted).

[Kai-Striega]

@mdhaber I'm not really sure how to do this, could you expand on it a little more?

[mdhaber]

By "grids", it sounds like the user is referring to intervals of integration, as defined by a and b, and the request is to vectorize quad w.r.t. a and b:

import numpy as np
import scipy.integrate

def new_quad(f, a, b):
    def g(a, b):
        return scipy.integrate.quad(f, a, b)
    h = np.vectorize(g)
    res = h(a, b)
    # np.asarray ensures that we don't return 0d arrays
    # tuple is in case we care about returning a tuple
    return  tuple(np.asarray(res))

def f(x):
    return x
a = 0
b = [1, 2, 3]
new_quad(f, a, b)  # (array([0.5, 2. , 4.5]), array([5.55111512e-15, 2.22044605e-14, 4.99600361e-14]))

[mdhaber]

I missed that quad_vec already exists. Vectorization of solve_ivp and nquad have been requested elsewhere, so I'll close this.
Update: oops. That doesn't handle vector a and b.