3

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?

user1118321
  • 25,567
  • 4
  • 55
  • 86
user2077948
  • 133
  • 5
  • 1
    I'm short on coffee this morning, but are you sure you've got the desired independent variables mapped to the limits of integration inside `adaptIntegrate` ? Generally if you get NaNs it's for a "good" reason, such as in fact having `f(x,y)` returning NaN for certain values of x or y. – Carl Witthoft Feb 16 '13 at 13:45
  • 2
    The domain of `asin` is `[-1,1]` for real values. `asin(sqrt((pi*area)/x))` will exceed this when `x – James Feb 16 '13 at 14:05
  • I'm not quite sure I understand what you mean by mapping the right variables to the limits of integration inside `adaptIntegrate`. I think the values of x and y are ok because the function (fun) produces numerical values when I just run `fun(238)`, for example. It's the integration which seems to get thrown and produce NaNs but I don't know how to see what's causing the problem within the code as it runs. As I said, I think it's the fact that it's a triangular shape to integrate over (as R will integrate when I leave out "asin" in `fun`) - I think I need to manipulate adaptIntegrate's working? – user2077948 Feb 16 '13 at 14:09
  • @James, I thought it might be a problem with the domains of `asin`, but that's why I created `newdistances`. Where L <- 25 distances <- L*seq(0.005, 100, 0.05) plot(distances) plot(log(distances)) roots <- sqrt((pi*area)/distances) plot(roots,type="l") It looks like I should be able to use values newdistances[64] onwards? (Sorry for the awful code formatting, not sure how to post code in a comment) – user2077948 Feb 16 '13 at 14:13
  • 1
    Actually, James, I see what you're saying - my newdistances isn't even incorporated into the integration in any way. I'll sort that out, thanks. – user2077948 Feb 16 '13 at 16:11

1 Answers1

2

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
agstudy
  • 119,832
  • 17
  • 199
  • 261
  • Thank you agstudy, this is very helpful. I appreciate your taking the time to write all this down and explain what's going wrong. – user2077948 Feb 16 '13 at 16:19
  • This has thrown up another question for me now (my math is *very* rusty!) - so the function can only work within the domain [pi*area,Inf] because the numerator has to be equal to or smaller than the denominator. This much, I understand (and thank you for simplifying the function so elegantly!). Now that the integration would be working within the correct domain, though, what does the integral really mean? Before, I could imagine the integration over a triangular area but now I don't see how the limits/integration relate(s) to d, x or r any more. – user2077948 Feb 16 '13 at 17:21
  • 1
    To imagine the integrals..plot the function see ?curve function in R...the inetgaral is just the area limited by upper and lower limits... – agstudy Feb 16 '13 at 17:29
  • Thank you very much. I can't tell you how helpful your responses have been. It's really helped me to think things through and understand a bit better. – user2077948 Feb 16 '13 at 20:13
  • @user2077948 glad that I can help. don't forget to accept th eanswer ( by checking it) if you find it helpful of course..:) – agstudy Feb 16 '13 at 20:47