scipy dblquad 在简单二重积分中提供错误的结果

Question

我正在尝试在Python中计算一个简单的双定积分：平方[0,1] x [0,1]中的函数Max(0, (4-12x) + (6-12y))。

我们可以用 Mathematica 来做并得到准确的结果：

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

通过简单的蒙特卡罗模拟，我可以确认这个结果。但是，使用

scipy.integrate.dblquad

我得到的值为 0.0005772072907971，错误为 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) )

同意该函数不可导，但它是连续的，我认为这个特定的积分没有理由有问题。

我想也许

dblquad

有一个“选项”参数，我可以在其中尝试不同的数值方法，但我没有发现类似的东西。

那么，我做错了什么？

Answer 1

尝试不同的数值方法

考虑到迭代

quad

 在 Windows 上遇到的麻烦，这就是我的建议。将其更改为显式两步过程后，您可以用另一种方法替换 

quad

 之一，

romberg

 对我来说似乎是最好的选择。 

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)

这会在我的计算机上打印 1.1574073959987758，希望您也能得到相同的结果。