Performing a double integral over a matrix of limits

Question

I have recently been learning how to perform double integrals in python. This is what I am using:

myint = dblquad(lambda xa,xb: np.exp(-(xa-xb)**2),-np.inf,x1,lambda x: -np.inf, lambda x: x2)

and for testing purposes I have chosen x1 and x2 to be say 5 and 10. This seems to work.

But in reality, my x1 = [1,2,3,4,5] and x2 = [5,6,7,8,9] and I want the double integral to be performed over every combination of x1 and x2 i.e. a matrix. I could do this with 2 for loops I guess, but I thought there might be a better way.

So my question is just - how do I perform a double integration over a matrix of limits.

Thank you.

---

edit:

I got the following warning:

UserWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze 
the integrand in order to determine the difficulties.  If the position of a 
local difficulty can be determined (singularity, discontinuity) one will 
probably gain from splitting up the interval and calling the integrator 
on the subranges.  Perhaps a special-purpose integrator should be used.

Does this mean that it doesn't converge? I don't really understand the message.

When I plot:

y = exp(-(x-5)^2)

for example, it just looks like a gaussian curve, so there is no problem integrating over that right? Is the problem because of the double integral?

Thank you.

---

edit:

Ah, I see. Thanks Raman Shah, I understand the problem now.

Answer 1

Using itertools you can create an iterator of limits to walk over. This essentially is a double loop, but for more extensible as you can have an arbitrary number of inputs with itertools.product and you don't store all the limits at once:

import numpy as np
from scipy.integrate import dblquad
import itertools

f = lambda xa,xb: np.exp(-(xa-xb)**2)
intg = lambda (x1,x2): dblquad(f,-np.inf,x1,
                               lambda x:-np.inf, 
                               lambda x:x2)


X1 = np.arange(1,6)
X2 = np.arange(5,10)
for limit in itertools.product(X1,X2):
    print limit, intg(limit)

If you need more speed, you can look into the multiprocessing module for parallel computation since each process is independent.

Answer 2

Why not use pythons zip function to feed exactly the values from each tuple that you want treated as inputs to your double integral, and then use map/apply to operate on those discrete pairs