【数値情報処理a】第5回の答え

$y=f(x)$ と $x$ 軸ではさまれた領域の面積をモンテカルロ法で求めるプログラムを作って提出しなさい。 ただし、$y$ が負の領域も正の面積として求める必要がある。

求める面積は下図の赤い部分である。

$f(x)=-x^{4}-x^{3}+10x^{2}+x-2$ の場合

図の赤い部分は $f(x)>y>0$ の領域、または $0>y>f(x)$ の領域である。

判別関数は $(x,y)$ が図の赤い部分に入れば 1 を返し、それ以外は 0 を返すようにする。

f <- function(x, y) {
  fx <- -x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2
  if (fx > y && y > 0) return(1)
  if (fx < y && y < 0) return(1)
  return(0)
}

図形を囲む長方形の $[x_{1},x_{2}]$ は $f(x)=0$ の解の最小値と最大値を使う。

x1 <- -3.63423216189401010467
x2 <- 2.71669197614524948747

図形を囲む長方形の $[y_{1},y_{2}]$ は $f(x)$ のグラフから大体の値を決めるとよい。

y1 <- -2.1
y2 <- 35

以上を踏まえてプログラムは次のようになる。

report5-1.R

# レポート
# モンテカルロ法
 
# 判別関数
f <- function(x, y) {
  fx <- -x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2
  if (fx > y && y > 0) return(1)
  if (fx < y && y < 0) return(1)
  return(0)
}
 
# 点の個数
n <- 1000
 
# 図形を囲む長方形の座標
x1 <- -3.63423216189401010467
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 35
 
# n個の一様乱数
x <- runif(n, x1, x2)
y <- runif(n, y1, y2)
 
m <- 0
for (i in 1 : n) {
  # 図形内部の点を判別関数でカウント
  m <- m + f(x[i], y[i])
}
 
# 面積
s <- m / n * (x2 - x1) * (y2 - y1)
print(s)

絶対値を使うと判別関数は次のようにできる。

f <- function(x, y) {
  if (y < abs(-x ^ 4 - p * x ^ 3 + 10 * x ^ 2 + x - 2) && y > 0) return(1)
  return(0)
}

$f(x)$ を絶対値にとると $y<0$ の領域は必要ないので $[y_{1},y_{2}]$ は次のようにしてよい。

y1 <- 0
y2 <- 35

$f(x)=-x^{4}-x^{3}+10x^{2}+x-2$ の場合

fx <- -x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2

x1 <- -3.63423216189401010467
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 35

$f(x)=-x^{4}-2x^{3}+10x^{2}+x-2$ の場合

  fx <- -x ^ 4 - 2 * x ^ 3 + 10 * x ^ 2 + x - 2

x1 <- -4.26427744793721785754
x2 <- 2.32569007375234804158
y1 <- -2.1
y2 <- 58.3

$f(x)=-x^{4}-3x^{3}+10x^{2}+x-2$ の場合

  fx <- -x ^ 4 - 3 * x ^ 3 + 10 * x ^ 2 + x - 2

x1 <- -4.95934144117426356502
x2 <- 2.00000000000000000000
y1 <- -2.1
y2 <- 96.1

$f(x)=-x^{4}-4x^{3}+10x^{2}+x-2$ の場合

  fx <- -x ^ 4 - 4 * x ^ 3 + 10 * x ^ 2 + x - 2

x1 <- -5.70992218440589560657
x2 <- 1.72955653934345376932
y1 <- -2.1
y2 <- 155.5

$f(x)=-x^{4}-5x^{3}+10x^{2}+x-2$ の場合

  fx <- -x ^ 4 - 5 * x ^ 3 + 10 * x ^ 2 + x - 2

x1 <- -6.50612662688904652697
x2 <- 1.50382273473682359644
y1 <- -2.1
y2 <- 245.8