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

Question

我知道 cubature 等 R 包能够对许多积分执行数值积分。但是，在我的情况下，积分的限制不是整数，而是其他变量的函数。

例如，假设我有

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

我想计算 x 范围内 y < x < 2 和 y 范围内 0 < y < 1 的二重积分 (dy dx)。我将如何在 R 中执行此操作（即这可能使用 pracma 包吗？）

如果这还没有实现，我会感到惊讶，但在那种情况下，link 算法会有所帮助。

Answer 1

如果你只是添加（或者更确切地说是乘以）一个将函数的值设置为零的项在任何 (x,y)-tuple where (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

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

n-dimensional integration in R with limits as functions

c

r

integration

numerical