3

私はプログラミングに非常に慣れておらず、本質的に試行錯誤によって学習してきましたが、アプローチ方法がわからない問題に到達しました。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 の積分中に条件付き引数の評価を実行できないように見えるため、積分したい正確な領域のみの値を取得できません。誰かアイデアやアドバイスはありますか？

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 つの項がありますが、問題を引き起こす可能性があるのは 1 つだけです。

    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 に答える