3

我对编程很陌生，基本上是通过反复试验来学习的，但是遇到了一个我不知道如何解决的问题。我需要在 R 中的三角形区域上进行双重积分。由于通常的积分函数似乎无法处理这个问题，我尝试使用cubature package(*edited - 完整代码见下文)。

更新/编辑： 我一直在努力解决这个问题，但仍然遇到同样的问题。我知道我必须确保值在关于 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进行积分时执行条件参数的评估，因此我仍然无法仅获得我想要积分的确切区域的值。有没有人有任何想法或建议？

4

1 回答 1

2

我不认为后面的代码adaptIntegrate会帮助你发生什么。您可以在控制台中输入adaptIntegrate，您将获得代码。它本质上是对 C 算法的调用。

  1. 为了了解它发生了什么，我认为您需要先了解您集成的内容。尽量简化你的功能，看看他的定义域。

     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
于 2013-02-16T16:09:57.837 回答