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

Question

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

我们可以用 Mathematica 完成并得到准确的结果：

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

通过简单的 Monte Carlo 模拟，我可以确认这个结果。但是，使用 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 有一个 'options' 参数，我可以在其中尝试不同的数值方法，但我没有找到类似的东西。

那么，我做错了什么？

Answer 1

try different numerical methods

考虑到迭代 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，希望你也能得到同样的结果。

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

scipy dblquad providing the wrong result in simple double integral

integration

scipy