講義（山本裕樹）

この文書は読取専用です。文書のソースを閲覧することは可能ですが、変更はできません。もし変更したい場合は管理者に連絡してください。
====== 【数値情報処理a】第5回の答え ======
===== 課題 =====

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

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

{{ report_montecarlo.png?nolink |求める面積}}

----
==== 解答例 ====

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

=== ①領域で判別する方法 ===

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

判別関数は $(x,y)$ が図の赤い部分に入れば 1 を返し、それ以外は 0 を返すようにする。
<code rsplus>
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)
}
</code>

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

<code rsplus>
x1 <- -3.63423216189401010467
x2 <- 2.71669197614524948747
</code>

図形を囲む長方形の $[y_{1},y_{2}]$ は $f(x)$ のグラフから大体の値を決めるとよい。
<code rsplus>
y1 <- -2.1
y2 <- 35
</code>

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

<file rsplus report5-1.R[enable_line_numbers="true"]>
# レポート
# モンテカルロ法

# 判別関数
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)
</file>

=== ②絶対値を使う方法 ===

絶対値を使うと判別関数は次のようにできる。
<code rsplus>
f <- function(x, y) {
  if (y < abs(-x ^ 4 - p * x ^ 3 + 10 * x ^ 2 + x - 2) && y > 0) return(1)
  return(0)
}
</code>

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

----
==== ①領域で判別する方法の $x_{1},x_{2}$ と $y_{1},y_{2}$ の範囲 ====

$f(x)=-x^{4}-x^{3}+10x^{2}+x-2$ の場合
<code rsplus>
fx <- -x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2
</code>

<code rsplus>
x1 <- -3.63423216189401010467
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 35
</code>


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

<code rsplus>
  fx <- -x ^ 4 - 2 * x ^ 3 + 10 * x ^ 2 + x - 2
</code>

<code rsplus>
x1 <- -4.26427744793721785754
x2 <- 2.32569007375234804158
y1 <- -2.1
y2 <- 58.3
</code>

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

<code rsplus>
  fx <- -x ^ 4 - 3 * x ^ 3 + 10 * x ^ 2 + x - 2
</code>

<code rsplus>
x1 <- -4.95934144117426356502
x2 <- 2.00000000000000000000
y1 <- -2.1
y2 <- 96.1
</code>

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

<code rsplus>
  fx <- -x ^ 4 - 4 * x ^ 3 + 10 * x ^ 2 + x - 2
</code>

<code rsplus>
x1 <- -5.70992218440589560657
x2 <- 1.72955653934345376932
y1 <- -2.1
y2 <- 155.5
</code>

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

<code rsplus>
  fx <- -x ^ 4 - 5 * x ^ 3 + 10 * x ^ 2 + x - 2
</code>

<code rsplus>
x1 <- -6.50612662688904652697
x2 <- 1.50382273473682359644
y1 <- -2.1
y2 <- 245.8
</code>