Double integral with variable boundaries in python Scipy + sympy (?)

Question

The full mathematical problem is here.

Briefly I want to integrate a function with a double integral. The inner integral has boundaries 20 and x-2, while the outer has boundaries 22 and 30.

I know that with Scipy I can compute the double integral with scipy.integrate.nquad. I would like to do something like this:

def f(x, y):
    return (x ** 2 + y ** 2)
res = sp.integrate.nquad(f, [[22, 30], [20, x-2]])

Is it possible? Maybe using also sympy?

Answer 1

I solved with sympy:

from sympy import *

x, y = symbols("x y")
f = (x ** 2 + y ** 2)
res = integrate(f, (y, 20, x-2), (x, 22, 30))

Basically sympy.integrate is able to deal with multiple integrations, also with variable boundaries.

Answer 2

If you need the numerical integration and sympy is not an option. Then you could try something like the following. For this example it seems quick, but I have a suspicion you may run into problems in general, see how well it does for your use case.Perhaps this possibly imperfect answer will prompt someone to submit something better.

I use the fact that we can do the integrations one after the other, integrating out the y first, to get a function of x, then integrating that.

from scipy.integrate import quad

def integrand(x, y): 
    return (x ** 2 + y ** 2)

def y_integral(x):
    # Note scipy will pass args as the second argument
    # we can fiddle it to work correctly, but by symmetry we don't need to here.
    return quad(integrand, 20, x-2, args=(x))[0]

We then use this y_integral function as the result function of the inner integral.

res = quad(y_integral, 22, 30) 
print res 

You could wrap this in a function if you use it regularly.