14

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?

Community
  • 1
  • 1
gc5
  • 9,468
  • 24
  • 90
  • 151

2 Answers2

20

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.

gc5
  • 9,468
  • 24
  • 90
  • 151
  • 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`.) – Thomas Ahle Mar 16 '22 at 00:21
6

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.

camz
  • 605
  • 3
  • 15