0

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?

zeycus
  • 860
  • 1
  • 8
  • 20
  • 1
    When I run this using scipy 1.1.0 and numpy 1.15.4, the output is `dblquad: 1.1574073869306833. Error: 0.0000000000031299` – Warren Weckesser Nov 12 '18 at 15:16
  • @Warren thank you, very interesting, that would be correct. I see I am using numpy 1.14.3, scipy 1.1.0. Python 3.6.5 – zeycus Nov 12 '18 at 16:00
  • Are you sure the code that you show in the question is exactly the same as the code that gave you `0.0005772072907971`? – Warren Weckesser Nov 12 '18 at 16:01
  • Yes, I just recopied from this screen to a new file, and executed it. Same 0.000577 result. – zeycus Nov 12 '18 at 16:36
  • 2
    Are you on Windows? How did you install scipy? You might be hitting this bug: https://github.com/scipy/scipy/issues/6882; see also https://stackoverflow.com/questions/27270044/big-integration-error-with-integrate-nquad – Warren Weckesser Nov 12 '18 at 17:30
  • Yes, I am. Installed everything through Anaconda, when I have time I think I should reinstall to the newest versions, see if the bug remains. Thanks! – zeycus Nov 13 '18 at 08:16

1 Answers1

1

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.