在R中三角形区域的双重积分

Question

在R中三角形区域的双重积分

3

我对编程非常陌生，基本上是通过试错学习的，但我遇到了一个我不知道如何解决的问题。我需要在 R 中对三角形区域进行双重积分。由于通常的 integrate 函数似乎无法处理这个问题，我尝试使用 cubature 包（*编辑-请参见下面的完整代码）。

更新/编辑：我已经更多地研究了这个问题，仍然遇到了同样的问题。我理解我必须确保相对于 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积分的同时进行条件参数while的评估，因此我仍然无法获得仅针对我想要积分的确切区域的值。是否有人有任何想法或建议？

- user2077948

1

今天早上我缺咖啡了，但你确定你已经将所需的独立变量映射到adaptIntegrate内的积分限制吗？通常如果出现NaN，这是一个“好”的原因，例如实际上f(x,y)对于某些x或y的值返回NaN。 - Carl Witthoft

2

asin 的定义域为实数 [-1,1]。当 x<pi*area ~= 94 时，asin(sqrt((pi*area)/x)) 将超出此范围，即对于积分中的所有值都是如此。当 asin 获得其定义域之外的值时，它会返回 NaN。显然，fun(235) 将给出一个有效的结果，但它远远超出了积分尝试计算的范围。 - James

我不太确定你所说的将正确的变量映射到adaptIntegrate内部的积分限制是什么意思。我认为x和y的值是正确的，因为当我运行例如fun(238)时，函数(fun)会产生数值。似乎是积分出了问题并产生了NaN，但我不知道如何在代码运行时查看导致问题的原因。正如我所说的，我认为这是由于要积分的形状是三角形（当我在fun中省略"asin"时，R将进行积分）- 我认为我需要操纵adaptIntegrate的工作方式？ - user2077948

@James，我认为这可能与asin的域有关，但这就是我创建newdistances的原因。其中L <- 25 distances <- Lseq(0.005, 100, 0.05) plot(distances) plot(log(distances)) roots <- sqrt((piarea)/distances) plot(roots,type="l")看起来我应该能够使用newdistances [64]及其以后的值？（很抱歉代码格式糟糕，不确定如何在注释中发布代码） - user2077948

1

其实，詹姆斯，我明白你的意思了——我的newdistances根本没有被整合到任何地方。我会解决这个问题的，谢谢。 - user2077948

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- agstudy · Accepted Answer

我认为adaptIntegrate背后的代码不会帮助你理解发生了什么。你可以在控制台中输入adaptIntegrate，然后你就能看到这段代码。实际上，它本质上是调用了一个C算法。

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只在[-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