Double integration over triangular area in R
I am very new to programming and have essentially studied the process and errors, but got to the problem, I don't know how to approach. I need to do double integration over a triangular area in R. Since the normal integral function doesn't seem to handle this, I tried using cubature package
(* edited - see below for complete code).
Update / Edit: I've worked on this more and I'm still facing the same problem. I understand that I have to ensure that the values ββare within the appropriate bounds with regard to the computation of asin. However, this still does not address the fundamental problem of the triangular region. It might be clearer if I write my complete 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 doesn't seem to be able to evaluate the conditional argument when integrating over x, so I still can't get values ββjust for the exact scope that I want to integrate. Anyone have any ideas or tips?
source to share
I don't think the code adaptIntegrate
will help you with what will happen. You can enter the console adaptIntegrate
and you get the code. This is essentially a call to Algorithm C.
-
To understand what is going on, I think you need first to understand what you are integrating. Try to simplify your function to see its 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)) }
There are two terms on X, but only one can create a problem.
asin(sqrt((pi*area)/X))
asin
defined only between [-1,1], sqrt
defined only for positive numbers.
So here fun is defined in between [pi*area,INF]
and you need to integrate in that domain.
eg:
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
source to share