Integrate a function that returns a matrix numpy

Question

I want to do something similar to this question.

Nm, Mm = np.meshgrid(range(3), range(2))
y = lambda x: x*Nm + x*Mm

Then, y returns a 3x2 matrix. I want to integrate y from a to b, for instance we can choose a=0 and b=1. This means that the i, j component of the integrated matrix should be int from 0 to 1 of (x * i + x * j) dx. If I take into account one of the answers:

>>> a = [sin, cos]
>>> vectorize(quad)(a, 0, pi)

Clearly a is a list of functions, but what I have is a function that returns an array, which is different. I get:

res = np.vectorize(integrate.quad)(y, 0, 1)

error: Supplied function does not return a valid float.

How can I fix this? Thanks for your help

Edit:

The desired result is

res = np.empty((3,2))
for i in range(3):
    for j in range(2):
        res[i, j] = quad(lambda x: x*i + x*j, 0, 1)[0]

So is the desired result the same as doing `quad` 6 times with different values of `i` and `j`? For reference show us how you'd do this with explicit loops. — hpaulj, Jun 03 '17 at 17:32

score 2 · Accepted Answer · answered Jun 03 '17 at 17:53

An alternative to passing the i,j via a lambda is to use the args parameter of quad:

In [39]: def foo(x,i,j):
    ...:     return x*i + x*j
In [40]: i,j=1,2
In [41]: integrate.quad(lambda x: x*i+x*j, 0,1)
Out[41]: (1.5, 1.6653345369377348e-14)
In [42]: integrate.quad(foo, 0, 1, args=(i,j))
Out[42]: (1.5, 1.6653345369377348e-14)

This is still evaluating the function separately for each i,j pair.

You may be able to bury the i,j iteration in a np.vectorize function, but that doesn't improve speed. I advocate it only if you need the broadcasting help.

A similar question comes up periodically. People want to apply quad (or other scipy functions) for a range of functions or parameters. But assuming you want to give quad full power over integration steps and such, there isn't a substitute for iteration.

score 1 · Answer 2 · answered Jun 03 '17 at 18:57

What you want can be achieved this way :

Nm, Mm = np.meshgrid(range(3), range(2))

def f(m,n):
    def g(x):
        return m*x+n*x
    return g

fv=np.vectorize(f)

u=fv(Mm,Nm)

np.vectorize(quad)(u,0,1)

But no numpy acceleration will arise here. it's just a commodity to mimic numpy loopless style.

Integrate a function that returns a matrix numpy

2 Answers2