我是编程的新手,基本上是通过反复试验来学习的,但遇到了一个我不知道如何解决的问题。我需要在 R 中的三角形区域上进行双重积分。由于通常的积分函数似乎无法处理这个问题,我尝试使用 cubature package
(*已编辑 - 完整代码见下文).
更新/编辑: 我一直在努力解决这个问题,但仍然遇到同样的问题。我知道我必须确保值在 asin 计算的适当范围内。然而,这仍然没有解决三角形区域的根本问题。如果我在下面发布我的完整代码,也许会更清楚:
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()
从警告消息来看,R 似乎无法执行条件参数的评估而 对 x 进行积分,因此我仍然无法仅获取我想要的确切区域的值整合过来。有没有人有任何想法或建议?
最佳答案
我认为 adaptIntegrate
背后的代码不会帮助您解决所发生的事情。您可以在控制台中输入 adaptIntegrate
,您将获得代码。它本质上是对 C 算法的调用。
为了理解它发生了什么,我认为你需要先理解你整合了什么。尽量简化你的功能,看看他的定义域。
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)) }
这里是 X 上的 2 个项,但只有一个会产生问题。
asin(sqrt((pi*area)/X))
asin
仅在[-1,1]之间定义,sqrt
仅针对正数定义。
所以这里fun是在[pi*area,INF]
之间定义的,你必须要集成到这个域中。
例如:
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
关于r - R中三角形区域的双重积分,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14909644/