Reputation: 108563
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.
Upvotes: 19
Views: 20993
Reputation: 10323
The cubature package does 2D (and N-D) integration using an adaptive algorithm. It should outperform simpler approaches for most integrands.
Upvotes: 18
Reputation: 121077
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.
Upvotes: 8