R 中的 n 维积分，极限为函数

Question

我知道 R 软件包（例如

cubature

）能够对许多积分执行数值积分。然而，在我的情况下，积分的限制不是整数，而是其他变量的函数。

例如，假设我有

f(x, y) = 4*(x-y)*(1-y)，其中 0 < y < 1, y < x < 2

我想计算 y < x < 2 for the range of x and 0 < y < 1 for the range of y. How would I do this in R (i.e. is this possible using the

pracma

包上的二重积分 (dy dx)？)

如果这还没有实现，我会感到惊讶，但在这种情况下，算法的链接会很有帮助。

Answer 1

如果您只需添加（或者乘以）一个项，将函数的值在任何 (x,y) 元组（其中 (y >= x)）处设置为零，这实际上非常简单

library(pracma)
fun <- function(x, y) 4*(x-y)*(1-y)*(y < x)
integral2(fun, 0, 2, 0, 1, reltol = 1e-10)
#=============
$Q
[1] 2.833363

$error
[1] 2.394377e-11

积分函数不会接受无限界，但由于

y

被限制为 (0-1) 并且

x

大于

y

，那么

x

的下界也必须为 0，这就是为什么我就是这样设定限制的。

Answer 2

您可以将域分割为三个三角形（制作一张图片）。通过在列中给出这三个三角形的顶点来定义它们：

T1 <- cbind(c(0,0), c(1,0), c(1,1))
T2 <- cbind(c(1,0), c(2,0), c(1,1))
T3 <- cbind(c(1,1), c(2,1), c(2,0))

然后你就可以使用

SimplicialCubature

包了。首先定义数组中三个三角形的并集：

Domain <- abind::abind(T1, T2, T3, along=3)

定义被积函数：

f <- function(xy) 4*(xy[1]-xy[2])*(1-xy[2])

然后：

library(SimplicialCubature)
> adaptIntegrateSimplex(f, Domain)
$integral
[1] 2.833333

$estAbsError
[1] 2.833333e-12

$functionEvaluations
[1] 96

积分的精确值为

17/6 = 2.833333....

。因此，我们得到的近似值比另一个答案中

pracma

给出的近似值更好（毫不奇怪：带有

(y<x)

项的被积函数并不平滑）。

---

对于这个例子，我们甚至可以使用

SimplicialCubature

得到更好的东西。被积函数是一个多项式：

f(x,y) = 4*x - 4*xy - 4*y + 4*y²

。

SimplicialCubature

包有多项式的精确方法。

> p <- definePoly(c(4,-4,-4,4), rbind(c(1,0),c(1,1),c(0,1),c(0,2)))
> printPoly(p)
Polynomial has 4 terms and 2 variables
4 * x[1]^1       (degree= 1 )
  -4 * x[1]^1 * x[2]^1       (degree= 2 )
  -4 * x[2]^1       (degree= 1 )
  + 4 * x[2]^2       (degree= 2 )
> integrateSimplexPolynomial(p, Domain)
$integral
[1] 2.833333

$functionEvaluations
[1] 9