rintegrationnumeric

n-dimensional integration in R with limits as functions


I'm aware of R packages such as cubature which are able to perform numerical integration for many integrals. However, in my situation the limits of integration are not integers, but functions of the other variables.

For instance, suppose I have

f(x, y) = 4*(x-y)*(1-y), where 0 < y < 1, y < x < 2

I want to compute a double integral (dy dx) over y < x < 2 for the range of x and 0 < y < 1 for the range of y. How would I do this in R (i.e. is this possible using the pracma package?)

I would be surprised if this isn't already implemented, but in that case, a link to an algorithm would be helpful.


Solution

  • It's actually quite easy if you just add (or rather multiply by) a term that sets the value of the function to zero at any (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
    

    The integral function will not accept an infinite bound (unlike the base R integrate function) , but since y is restricted to (0-1) and x is greater than y, then x's lower bound must also be 0, which was why I put the limits in that way.