4

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.

nathanesau
  • 1,681
  • 16
  • 27

2 Answers2

4

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 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.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • The condition should be ``x >= y`` but thanks for this answer – nathanesau Sep 19 '15 at 14:56
  • I actually posted the answer where I used `*(y < x)`. Then I looked at it and decided it was flipped and changed it , but now that I review the logic I do agree you are right and will revert. – IRTFM Sep 19 '15 at 19:45
0

You can split the domain into three triangles (make a picture). Define these three triangles by giving their vertices in columns:

T1 <- cbind(c(0,0), c(1,0), c(1,1))
T2 <- cbind(c(1,0), c(2,0), c(1,1))
T3 <- cbind(c(1,1), c(2,1), c(2,0))

Then you can use the SimplicialCubature package. First define the union of the three triangles in an array:

Domain <- abind::abind(T1, T2, T3, along=3)

Define the integrand:

f <- function(xy) 4*(xy[1]-xy[2])*(1-xy[2])

Then:

library(SimplicialCubature)
> adaptIntegrateSimplex(f, Domain)
$integral
[1] 2.833333

$estAbsError
[1] 2.833333e-12

$functionEvaluations
[1] 96

The exact value of the integral is 17/6 = 2.833333..... Thus the approximation we get is better than the one given by pracma in the other answer (not surprising: the integrand with the term (y<x) is not smooth).

For this example, we can even get something better with SimplicialCubature. The integrand is a polynomial: f(x,y) = 4*x - 4*xy - 4*y + 4*y². The SimplicialCubature package has an exact method for polynomials.

> p <- definePoly(c(4,-4,-4,4), rbind(c(1,0),c(1,1),c(0,1),c(0,2)))
> printPoly(p)
Polynomial has 4 terms and 2 variables
4 * x[1]^1       (degree= 1 )
  -4 * x[1]^1 * x[2]^1       (degree= 2 )
  -4 * x[2]^1       (degree= 1 )
  + 4 * x[2]^2       (degree= 2 )
> integrateSimplexPolynomial(p, Domain)
$integral
[1] 2.833333

$functionEvaluations
[1] 9
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225