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?

Susannahsusanne

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.

Susannahsusanne

Note that variables have to appear in the upper part of the range, not the lower. Something like res = integrate(f, (y, 20, 22), (x, y, 30)) doesn't work. (It returns a function in y.) – Happy

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.

Hexameter