Integrate a function that takes as input a scalar value and outputs a matrix

Question

I have a function that takes as an input a scalar value, and outputs a matrix (following is an MRE with size 2x2, my actual function has problem specific matrix dimensions, like say 100x100). I'd like to integrate said function over a range.

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

quad(func, 0, np.pi)

This gives me the following error :

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/scipy/integrate/_quadpack_py.py", line 465, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/scipy/integrate/_quadpack_py.py", line 577, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: only size-1 arrays can be converted to Python scalars

I wanted to avoid looping over each individual function, so I vectorized the integration operation, in the hopes that it would accommodate this sort of a function :

grid_quad = np.vectorize(quad)
grid_quad(func, 0, np.pi)

I still get the same error :

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/numpy/lib/function_base.py", line 2329, in __call__
    return self._vectorize_call(func=func, args=vargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/numpy/lib/function_base.py", line 2407, in _vectorize_call
    ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/numpy/lib/function_base.py", line 2367, in _get_ufunc_and_otypes
    outputs = func(*inputs)
              ^^^^^^^^^^^^^
  File "/scipy/integrate/_quadpack_py.py", line 465, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/scipy/integrate/_quadpack_py.py", line 577, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: only size-1 arrays can be converted to Python scalars

When I explicitly wrote out my integrand as an array of arrays, this method worked

efunc = [[np.cos, np.sin],[np.sin, np.cos]]
grid_quad(efunc, 0, np.pi)

and gave the expected output

(array([[4.92255263e-17, 2.00000000e+00],
       [2.00000000e+00, 4.92255263e-17]]), array([[2.21022394e-14, 2.22044605e-14],
       [2.22044605e-14, 2.21022394e-14]]))

which is a matrix that contains the values obtained by integrating the "function" associated with each matrix position.

In pen and paper terms, I have a matrix whose entries are functions that depend on a variable. I would like to integrate each one of these functions over an interval, and receive as output a matrix of the definite integral values.

Is there a "(num)pythonic" way to achieve this?

Answer 1

Try scipy.integrate.quad_vec.

import numpy as np
from scipy.integrate import quad_vec

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

quad_vec(func, 0, np.pi)
# (array([[2.22044605e-16, 2.00000000e+00],
#         [2.00000000e+00, 2.22044605e-16]]),
#  1.3312465980656846e-13)

The first output is your result, the second is the error estimate.