scipy dblquad 在简单的双积分中提供了错误的结果答案

【问题标题】：scipy dblquad providing the wrong result in simple double integralscipy dblquad 在简单的双积分中提供了错误的结果
【发布时间】：2019-04-15 07:55:03
【问题描述】：

我正在尝试在 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 有一个“选项”参数，我可以在其中尝试不同的数值方法，但我没有发现类似的东西。

那么，我做错了什么？

【问题讨论】：

当我使用 scipy 1.1.0 和 numpy 1.15.4 运行它时，输出为 dblquad: 1.1574073869306833. Error: 0.0000000000031299
@Warren 谢谢你，非常有趣，那是正确的。我看到我正在使用 numpy 1.14.3，scipy 1.1.0。 Python 3.6.5
您确定您在问题中显示的代码与提供给您0.0005772072907971 的代码完全相同吗？
是的，我刚刚从这个屏幕复制到一个新文件，并执行了它。相同的 0.000577 结果。
您在 Windows 上吗？你是如何安装 scipy 的？你可能遇到了这个错误：github.com/scipy/scipy/issues/6882;另见stackoverflow.com/questions/27270044/…

标签： scipy integration

【解决方案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，希望你也能得到同样的结果。

【讨论】：

这确实解决了，谢谢！我不知道romberg 函数。