使用 Python 进行双重集成答案

【问题标题】：Double integration using Python使用 Python 进行双重集成
【发布时间】：2020-12-01 18:06:42
【问题描述】：

我有以下积分（我不知道为什么屏幕截图这么低，所以你需要滚动一下）。我尝试了以下代码：

from scipy import pi, sqrt
from scipy.integrate import dblquad
func = lambda x,y: ((1-sqrt(x**2+y**2))/(sqrt(x**2+y**2)))*(x**2)
x1,x2 = -1, 1
y1,y2 = lambda x: -sqrt(1-x**2), lambda x: sqrt(1-x**2)
print (dblquad(func, x1, x2, y1, y2))

但我不断得到

The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  **opt)

我得到的结果似乎不正确：

(nan, 1.0867783627477699e-07)

【问题讨论】：

我怀疑问题来自于 SciPy 在 x = 0、y = 0 评估您的函数，它会给出一个 nan。尝试重写函数以在这种特殊情况下返回 0，例如类似func = lambda x, y: 0.0 if x == y == 0 else ((1-sqrt(x**2+y**2))/(sqrt(x**2+y**2)))*(x**2).
@MarkDickinson 谢谢，它有效！

标签： python math scipy

【解决方案1】：

积分的被积函数未在点 (0, 0) 处定义。从数学上讲，你最终将0 除以0。在计算上，func(0.0, 0.0) 会给出一个 NaN，所以如果 dblquad 恰好在此时计算你的被积函数，你可以期望得到一个 NaN 作为结果。

然而，被积函数在(0, 0)(0, 0) 附近是连续的，并且在(x, y) 从任何方向接近(0, 0) 时具有0 的限制，因此扩展定义的定义是很自然的通过指定其在(0, 0) 处的值为0 来对整个单元盘进行积分。

因此，要修复集成，您可以将 func 重写为特例点 (0, 0)，例如：

def func(x, y):
    if x == y == 0.0:
        return 0.0
    else:
        return (1-sqrt(x**2+y**2))/(sqrt(x**2+y**2)))*(x**2)

在我的机器上，如果我对func 进行替换，然后完全按照编写的方式执行其余代码，我会得到以下结果：

>>> print (dblquad(func, x1, x2, y1, y2))
(0.261799387785821, 1.4745421123040386e-08)

而且我相信您知道，通过分析评估积分以发现值应该是 π/12 很简单。那么让我们比较一下我们得到的值和期望值：

>>> pi / 12.
0.2617993877991494

所以看起来 SciPy 给出的结果在所需的误差范围内是准确的。

【讨论】：