在R中快速计算二重积分

Question

我正在寻找一个比

<scipy.integrate.nquad>

更快的二重积分解决方案。

integrate(function(y) { 
   sapply(y, function(y) {
     integrate(function(x) myfun(x,y), llim, ulim)$value
   })
 }, llim, ulim)

举个例子

myfun <- function(x,y) cos(x+y)
llim <- -0.5
ulim <- 0.5

我发现了一篇参考了名为quad2d的FORTRAN程序的古老论文，但是除了MATLAB的一些帮助页面外，我找不到其他的资源。因此，我正在寻找一个可以快速执行双重积分（即无需使用sapply循环）的C或FORTRAN库，并且可以从R中调用。只要它们都与GPL兼容，所有想法都非常值得欢迎。

如果解决方案涉及调用已经随R一起发布的库中的其他函数，我也很乐意听取这些方案。

- Joris Meys

---

你考虑过使用pracma::dblquad、pracma::simpson2d以及cubature和R2Cuba包中的函数吗？ - Joshua Ulrich

Answer 1

cubature软件包使用自适应算法进行 2D（和 N-D）积分。 对于大多数被积函数，该软件包应该优于简单的方法。

- Jeffrey Sax

---

谢谢指引，我会在这个周末运行测试，看哪一个让我最满意。 - Joris Meys

供今后参考，使用“cubature”方法的代码如下：adaptIntegrate(function(x) cos(x[1] + x[2]), c(llim, llim), c(ulim, ulim))。 - jbaums

嗨，我也在使用integrate进行双重积分时遇到了困难，我发现这很慢，但我不知道如何使用adaptIntegrate指定积分限制。我想对第一个积分在35和u之间进行积分，并对第二个积分在35和95之间进行积分，其中u在35和95之间移动。你会怎么写呢？谢谢。 - Mbr Mbr

Answer 2

Joshua提到的pracma软件包包含quad2d的一个版本。

quad2d(myfun, llim, ulim, llim, ulim)

这将给出与您的循环相同的答案，但要在数值容差范围内，使用示例函数。

默认情况下，对于您的示例函数，quad2d 比循环慢。如果您将 n 降低，可以加快速度，但我猜这取决于您的函数平滑程度以及您愿意为速度牺牲多少准确性。

- Richie Cotton

---

谢谢你提供的示例。我会在真实的函数中进行测试，因为它们有点复杂。准确性对我来说最重要，但速度确实是一个问题。我会随时向你更新进展情况。 - Joris Meys