来自 R 包 cubature 的 adaptIntegrate 得到一个巨大的积分错误

Question

我对各种 2D 求积工具进行了一些测试，cubature 中的 adaptIntegrate 一直是最精确的工具之一，当它工作时。

麻烦的是，在某些情况下它完全出错了，但不是小数点，真的超出了规模，我不明白为什么。

我试图在矩形（实际上甚至是正方形）二维域上集成的功能是合乎逻辑的，形式如下：

条件_1 & (条件_2 | 条件_3 | ...)

这对我尝试过的任何工具都没有问题（pracma::integral2、pracma::quad2d、pracma::simpson2d...）。
adaptIntegrate 然而，在大多数情况下提供最佳结果的同时，偶尔会完全失败。

例子：

require(cubature)

# Intersection of a circle of radius 4 centred in (0,0) with a circle of radius 1 centred in (0,0)

adaptIntegrate(function(x) ( (x[1]^2+x[2]^2 <= 16) & ((x[1]-0)^2+(x[2]-0)^2 <= 1) ), c(-4, -4), c(4, 4), absError = 1.e-2 )

#$integral
#[1] 3.141522
#
#$error
#[1] 0.009982141
#
#$functionEvaluations
#[1] 24089
#
#$returnCode
#[1] 0

正确：积分应该是 pi*1^2，大约是 3.142。
现在将小圆圈的中心移动到 (1,0)。它仍然完全包含在较大的圆圈中，因此交叉点仍然具有相同的区域。

# Intersection of a circle of radius 4 centred in (0,0) with a circle of radius 1 centred in (1,0)

adaptIntegrate(function(x) ( (x[1]^2+x[2]^2 <= 16) & ((x[1]-1)^2+(x[2]-0)^2 <= 1) ), c(-4, -4), c(4, 4), absError = 1.e-2 )

#$integral
#[1] 0
#
#$error
#[1] 0
#
#$functionEvaluations
#[1] 51
#
#$returnCode
#[1] 0

我无法弄清楚为什么会发生这种情况。如果我对小圆圈的 x 使用 1.1 而不是 1，它会回到一个非常令人满意的估计值。

如果我做错了什么，或者这实际上是一个错误，有没有人有任何想法？

谢谢！

PS
2 个注释：

1. 是的，我知道这些例子很简单，我可以简单地使用小圆圈的面积。我需要解决的真正问题不那么微不足道（例如，小圆圈可以部分重叠大圆圈的周边，并且会有很多小圆圈，不仅是一个，它们之间也会重叠）。如果我不能依靠它来处理一个简单的案例，我就不能依靠它来处理更复杂的案例，可以吗？
2. 是的，我看到了一些与adaptIntegrate 不一致的帖子，例如this one。
   我在那里看到了这个答案：

Cucture C 库给出的结果与我在上述问题中报告的结果相同，而且这不太可能是错误；相反，在某些情况下，h 自适应 cubature 例程（R 包是其接口）不如 Cubature 的 p 自适应例程准确，后者将适当区域的采样点数加倍。

这有什么帮助？当工具实际上给出错误答案时，我说 “这不是错误” 对我来说似乎是不正确的。此外，在我的情况下，pcubature 一直不太好，并且给出了相同类型的错误，所以......

Answer 1

我无法解释为什么adaptIntegrate无法返回合理的结果。你可以写信给cubature的维护者或者这个软件的原作者。

但是，如果您将积分域限制为围绕较小/最小圆的矩形，您可以获得合理的结果！

最好将函数定义为返回数字 价值观。这也更容易写成几个圈子。

fn10 <- function(xy) {
    x <- xy[1]; y <- xy[2]
    (x^2+y^2 <= 16) * ((x-1)^2 + y^2 <= 1)
}

fn40 <- function(xy) {
    x <- xy[1]; y <- xy[2]
    (x^2+y^2 <= 16) * ((x-4)^2 + y^2 <= 1)
}

library(cubature)

adaptIntegrate(fn10, lowerLimit = c(0, -1), upperLimit = c(2, 1))
## $integral
## [1] 3.141828

adaptIntegrate(fn40, lowerLimit = c(3, -1), upperLimit = c(5, 1))
## $integral
## [1] 1.490316

我产生了兴趣并尝试通过不同的集成例程为f40 区域找到更好的价值。

library(pracma)

fun40 <- function(x, y)
    (x^2+y^2 <= 16) * ((x-4)^2 + y^2 <= 1)

quad2d(fun40, 3, 5, -1, 1, n = 512)                ## 200 ms
## [1] 1.489193

integral2(fun40,  3, 5, -1, 1)                     ##   3 ms
## [1] 1.494751

simpson2d(fun40, 3, 5, -1, 1, nx = 512, ny = 512)  ##  16 ms
## [1] 1.487213

那么这个积分的真实值是多少？ 两个圆相交于x0 = 31/8，因此一维 在区间[3, 4] 上积分的（矢量化）函数是：

x0 <- 31/8.0
fun <- function(x) {
    ifelse(x <= x0, sqrt(1 - (x-4)^2), sqrt(16 - x^2))
}

2 * quadgk(fun, 3, 4)
## [1] 1.487332

我们将积分加倍，因为我们只在 x 轴上方积分。 这应该精确到大约 0.00025 的误差。我们看到 简单的 Simpson 2D 求积非常接近。

注意：我们也可以利用三角公式来计算面积，即两个圆的面积，因此

# circle at (0.0)               circle at (4,0)
th1 = 2*acos(x0/4);             th2 = 2*acos((4-x0)/1)
A1 = 0.5*(th1 - sin(th1))*16;   A2 = 0.5*(th2 - sin(th2))*1

A1 + A2
## [1] 1.487332

这个值与之前在一维积分的帮助下找到的值一致。