# [Solved] Integrating functions that return an array in Python

I have a lot of data to integrate over and would like to find a way of doing it all with just matrices, and would be willing to compromise on accuracy for a performance boost. What I have in mind is something like this:

``````import numpy
import scipy

a = np.array([1,2,3])

def func(x):
return x**2 + x

def func2(x):
global a
return a*x

def integrand(x):
return func(x)*func2(x)

``````

So I am trying to integrate each element in the array that comes out of `integrand`.

I’m aware that there is a possibility of using `numpy.vectorize()` like this:

``````integrated = numpy.vectorize(scipy.integrate.quad)(integrand, 0, 1)
``````

but I can’t get that working. Is there a way to do this in python?

Solution

Well now that I learnt a bit more python I can answer this question if anyone happens to stable upon it and has the same question. The way to do it is to write the functions as though they are going to take scalar values, and not vectors as inputs. So follow from my code above, what we would have is something like

``````import numpy as np

a = np.array([1, 2, 3]) # arbitrary array, can be any size

def func(x):
return x**2 + x

def func2(x, a):
return a*x

def integrand(x, a):
return func(x)*func2(x, a)

def integrated(a):
integrated, tmp = scipy.integrate.quad(integrand, 0, 1, args = (a))
return integrated

def vectorizeInt():
global a
integrateArray = []
for i in range(len(a)):
integrate = integrated(a[i])
integrateArray.append(integrate)
return integrateArray
``````

Not that the variable which you are integrating over must be the first input to the function. This is required for scipy.integrate.quad. If you are integrating over a method, it is the second argument after the typical `self` (i.e. `x` is integrated in `def integrand(self, x, a):`). Also the `args = (a)` is necessary to tell `quad` the value of `a` in the function `integrand`. If `integrand` has many arguments, say `def integrand(x, a, b, c, d):` you simply put the arguments in order in `args`. So that would be `args = (a, b, c, d)`.

## Solution #1:

`vectorize` won’t provide help with improving the performance of code that uses `quad`. To use `quad`, you’ll have to call it separately for each component of the value returned by `integrate`.

For a vectorized but less accurate approximation, you can use `numpy.trapz` or `scipy.integrate.simps`.

Your function definition (at least the one shown in the question) is implemented using numpy functions that all support broadcasting, so given a grid of `x` values on [0, 1], you can do this:

``````In [270]: x = np.linspace(0.0, 1.0, 9).reshape(-1,1)

In [271]: x
Out[271]:
array([[ 0.   ],
[ 0.125],
[ 0.25 ],
[ 0.375],
[ 0.5  ],
[ 0.625],
[ 0.75 ],
[ 0.875],
[ 1.   ]])

In [272]: integrand(x)
Out[272]:
array([[ 0.        ,  0.        ,  0.        ],
[ 0.01757812,  0.03515625,  0.05273438],
[ 0.078125  ,  0.15625   ,  0.234375  ],
[ 0.19335938,  0.38671875,  0.58007812],
[ 0.375     ,  0.75      ,  1.125     ],
[ 0.63476562,  1.26953125,  1.90429688],
[ 0.984375  ,  1.96875   ,  2.953125  ],
[ 1.43554688,  2.87109375,  4.30664062],
[ 2.        ,  4.        ,  6.        ]])
``````

That is, by making `x` an array with shape (n, 1), the value returned by `integrand(x)` has shape `(n, 3)`. There is one column for each value in `a`.

You can pass that value to `numpy.trapz()` or `scipy.integrate.simps()`, using `axis=0`, to get the three approximations of the integrals. You’ll probably want a finer grid:

``````In [292]: x = np.linspace(0.0, 1.0, 101).reshape(-1,1)

In [293]: np.trapz(integrand(x), x, axis=0)
Out[293]: array([ 0.583375,  1.16675 ,  1.750125])

In [294]: simps(integrand(x), x, axis=0)
Out[294]: array([ 0.58333333,  1.16666667,  1.75      ])
``````

Compare that to repeated calls to `quad`:

``````In [296]: np.array([quad(lambda t: integrand(t)[k], 0, 1)[0] for k in range(len(a))])
Out[296]: array([ 0.58333333,  1.16666667,  1.75      ])
``````

Your function `integrate` (which I assume is just an example) is a cubic polynomial, for which Simpson’s rule gives the exact result. In general, don’t expect `simps` to give such an accurate answer.

## Solution #2:

quadpy (a project of mine) is fully vectorized. Install with

``````pip install quadpy
``````

and then do

``````import numpy

def integrand(x):
return [numpy.sin(x), numpy.exp(x)]  # ,...

``````[0.45969769 1.71828183]