講義（山本裕樹）

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

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

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

{{ 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.63423216189423481381
x2 <- 2.71669197614524948747
</code>

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

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

<file rsplus report4.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.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)
</file>

----

$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.63423216189423481381
x2 <- 2.71669197614524948747
y1 <- -2.1
y2 <- 34.9
</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.26427744793721430483
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.00000000000002620126
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.70992218440590715289
x2 <- 1.72955653934345376932
y1 <- -2.1
y2 <- 155.4
</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.50612662688896659091
x2 <- 1.50382273473682359644
y1 <- -2.1
y2 <- 245.8
</code>


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

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

<code rsplus>
x1 <- -7.33896392537833719416
x2 <- 1.31293708037980039194
y1 <- -2.1
y2 <- 378.7
</code>