Search code examples
pythonscipynumerical-integration

scipy dblquad providing the wrong result in simple double integral

I am trying to calculate a straightforward doble definite integral in Python: function Max(0, (4-12x) + (6-12y)) in the square [0,1] x [0,1].

We can do it with Mathematica and get the exact result:

Integrate[Max[0, 4-12*u1 + 6-12*u2], {u1, 0, 1}, {u2, 0,1}] = 125/108.

With a simple Monte Carlo simulation I can confirm this result. However, using scipy.integrate.dblquad I am getting a value of 0.0005772072907971, with error 0.0000000000031299

from scipy.integrate import dblquad

def integ(u1, u2):
    return max(0, (4 - 12*u1) + (6 - 12*u2))
sol_int, err = dblquad(integ, 0, 1, lambda _:0, lambda _:1, epsabs=1E-12, epsrel=1E-12)
print("dblquad: %0.16f. Error: %0.16f" % (sol_int, err) )

Agreed that the function is not derivable, but it is continuous, I see no reason for this particular integral to be problematic.

I thought maybe dblquad has an 'options' argument where I can try different numerical methods, but I found nothing like that.

So, what am I doing wrong?

Solution

  • try different numerical methods

    That's what I would suggest, given the trouble that iterated quad has on Windows. After changing it to an explicit two-step process, you can replace one of quad with another method, romberg seems the best alternative to me. 

    from scipy.integrate import quad, romberg

def integ(u1, u2):
    return max(0, (4 - 12*u1) + (6 - 12*u2))

sol_int = romberg(lambda u1: quad(lambda u2: integ(u1, u2), 0, 1)[0], 0, 1)    
print("romberg-quad: %0.16f " % sol_int)

    This prints 1.1574073959987758 on my computer, and hopefully you will get the same.