19

I'm looking for a solution for a double integral that is faster than 

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

with eg

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

I found an old paper that referred to a FORTRAN program called quad2d, but I couldn't find anything else but some help pages for matlab for the rest. So I'm looking for a C or FORTRAN library that can do double integrals quick (i.e. without the sapply loop), and that can be called from R. All ideas are very much appreciated, as long as they're all GPL compatible.

If the solution involves calling other functions from the libraries that are already shipped with R, I'd love to hear from them as well.

Joris Meys
  • 106,551
  • 31
  • 221
  • 263

2 Answers2

18

The cubature package does 2D (and N-D) integration using an adaptive algorithm. It should outperform simpler approaches for most integrands.

Jeffrey Sax
  • 10,253
  • 3
  • 29
  • 40
  • Thx for the pointer, I'm running the tests this weekend and see which one is making me most happy. – Joris Meys Jan 18 '12 at 21:27
  • 6
    For future reference, the `cubature` approach would be: `adaptIntegrate(function(x) cos(x[1] + x[2]), c(llim, llim), c(ulim, ulim))`. – jbaums Jul 15 '14 at 01:06
  • Hi, I am also struggling with a double integral with integrate which I find quite slow but I don't know understand how you specify the limits of integration with adaptIntegrate I want to integrate between 35 and u for the first integral and between 35 and 95 for the second, with u moving between 35 and 95. How would you write it ? Thanks – Mbr Mbr May 06 '22 at 14:32
8

The pracma package that Joshua pointed out contains a version of quad2d.

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

This gives the same answer, within numerical tolerance, as your loop, using the example function.

By default, with your example function, quad2d is slower than the loop. If you drop n down, you can make it faster, but I guess it depends upon how smooth your function is, and how much accuracy you are willing to sacrifice for speed.

Richie Cotton
  • 118,240
  • 47
  • 247
  • 360
  • 1
    Thx for the example as well. I'll test it on the real functions, as they're a bit more complex. I need accuracy the most, but speed really is an issue. I'll keep you updated. – Joris Meys Jan 18 '12 at 21:29