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

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

• 求める面積は下図の赤い部分である。
• 第2回課題の4つの解のうち、最小値を $x_{1}$、最大値を $x_{2}$ として使うこと。

$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.63423216189423481381
x2 <- 2.71669197614524948747

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

y1 <- -2.1
y2 <- 34.9

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

report4.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.63423216189423481381
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 34.9
 
# 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(x)=-x^{4}-x^{3}+10x^{2}+x-2$ の場合

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

x1 <- -3.63423216189423481381
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 34.9

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

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

x1 <- -4.26427744793721430483
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.00000000000002620126
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.70992218440590715289
x2 <- 1.72955653934345376932
y1 <- -2.1
y2 <- 155.4

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

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

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

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

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

x1 <- -7.33896392537833719416
x2 <- 1.31293708037980039194
y1 <- -2.1
y2 <- 378.7