Double integration over a triangular area in R

Question

I am very new to programming and have been essentially learning by trial and error, but have reached a problem I do not know how to approach. I need to do a double integration over a triangular area in R. As the usual integrate function doesn't seem able to handle this, I tried using cubature package (*edited - see below for the full code).

Update/Edit: I've been working on this more and am still coming up against the same issue. I understand that I have to ensure that values are within the appropriate bounds with respect to the asin calculation. However, this still isn't getting around the fundamental problem of the triangular area. Perhaps it will be clearer if I post my full code below:

L <- 25
n <- -4
area <- 30
distances <- L*seq(0.005, 100, 0.05)

cond <- area*pi
d <- 5

fun <- function(x=1,r=0)
{
  if (x<cond) {
    return(0)
  } else {
    return((-1)*((n+2)/(2*pi*(L^2)))*(1+((x/L)^2))^(n/2)*(1/pi)*(1/pi)*acos(d/x))*asin(sqrt((pi*area)/d+r))
  }
}
fun(5)
fun(300)

library(cubature)
integrationone <- function()
{
  integrand <- adaptIntegrate(fun, lowerLimit=c(d,0), upperLimit=c(80,80))
  return(integrand$integral)
}
integrationone()
warnings()

From looking at the warning messages, R seems unable to carry out the evaluation of the conditional argument while integrating over x, so I still can't get values for only the exact area I want to integrate over. Does anyone have any ideas or advice?

Answer 1

I don't think that the code behind adaptIntegrate will help you what's happen. You can type in a console adaptIntegrate and you will get the code. It is essentially a call to a C algorithm.

1. In order to understand what it is happen , I think you need before to understand what you integrate. Try to simplify your function, to see his definition domain.

    INV_PI <- 1/pi
    fun <- function(X){  
       scale <- -1*((n+2)/(2*pi*(L^2)))*INV_PI^2 *acos(d/(d+r))
       res <- scale*asin(sqrt((pi*area)/X))* (1+((X/L)^2))^(n/2)
       sqrt(prod(res))
    }
   

   Here the 2 terms on X , but only one can produce problem.

   asin(sqrt((pi*area)/X)) 
   

asin is defined only between[-1,1], sqrt is defined only for positive numbers.

So here fun is defined between [pi*area,INF], and you have to integrate in this domain.

for example :

low.Lim <- pi*area
doubleintegration <- function()  
{  
  integrand <- adaptIntegrate(fun, lowerLimit=c(low.Lim,low.Lim), 
                                   upperLimit=c(200*low.Lim,200*low.Lim))  
  return(integrand$integral)  
}  

doubleintegration()
[1] 0.1331089