R语言含无穷限与函数边界的二重积分求解方法问询

你的需求是x从0到正无穷、内层y的上限为x的二重积分，可以将其拆解为累次积分，利用R基础包的integrate()函数（支持无穷积分限）嵌套实现，无需额外安装工具包。

具体实现步骤 #

1. 定义区分外层变量x和内层变量y的被积函数；
2. 编写内层积分函数：对给定的x，计算y从0到x的积分；
3. 对内层积分结果，执行x从0到正无穷的外层积分。

针对你的示例代码，完整实现如下： #

# 定义被积函数，x为外层变量，y为内层变量
f <- function(x, y) {
  exp(-(1/2)*x - (1/3)*y)
}

# 内层积分：对每个x，计算y从0到x的积分
inner_integral <- function(x) {
  # 提取integrate返回结果中的积分值
  integrate(function(y) f(x, y), lower = 0, upper = x)$value
}

# 外层积分：x从0到正无穷
result <- integrate(inner_integral, lower = 0, upper = Inf)

# 查看最终结果
print(result)

结果验证（可选） #

针对该被积函数，可通过解析解验证代码正确性：
二重积分$\int_{0}^{{\infty}\int_{0}}{x} e^{-\frac{x}{2}-\frac{y}{3}} dydx$的解析解为$\frac{12}{5}=2.4$，运行代码后会得到与该值近似的结果。

备选方案：使用cubature包的自适应积分 #

若偏好多变量积分函数，可通过cubature包的hcubature()函数，结合区域指示函数实现函数边界的积分：

install.packages("cubature")
library(cubature)

f <- function(x) {
  # x[1]对应x，x[2]对应y；仅当0 ≤ y ≤ x时计算被积函数，否则返回0
  if (x[2] >= 0 && x[2] <= x[1]) {
    exp(-(1/2)*x[1] - (1/3)*x[2])
  } else {
    0
  }
}

# 积分区域设为x从0到Inf、y从0到Inf，通过指示函数限制y≤x
result <- hcubature(f, lowerLimit = c(0, 0), upperLimit = c(Inf, Inf))
print(result$integral)

内容的提问来源于stack exchange，提问作者shinramen