講義（山本裕樹）

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

$y=f(x)$ と $x$ 軸ではさまれた領域の面積をシンプソンの公式で数値積分するプログラムを作って提出しなさい。
ただし、$y$ が負の領域も正の面積として求める必要がある。

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

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

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

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

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

関数 $f(x)$ は絶対値 ''abs()'' を使って定義する。

積分区間には $f(x)=0$ の解のうち、端の二つを使う。
<code>
-3.63423216189401010467
 2.71669197614524948747
</code>

<file rsplus report4-1.R[enable_line_numbers="true"]>
# レポート
# シンプソンの公式
# 絶対値を使う方法

# 被積分関数の定義
f <- function(x) {
  return(abs(-x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2))
}

# 区間数
n <- 1000

# 区間
a <- -3.63423216189401010467
b <-  2.71669197614524948747

# 幅
h <- (b - a) / n

s1 <- f(a + h * 0.5)
s2 <- 0
for (i in 1 : (n - 1)) {
  x <- a + i * h
  s1 <- s1 + f(x + h * 0.5)
  s2 <- s2 + f(x)
}
s <- (f(a) + 4 * s1 + 2 * s2 + f(b)) * h / 6
print(sprintf("n = %d", n), quote = FALSE)
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("S = %f", s), quote = FALSE)
</file>

=== ②積分範囲を分ける方法 ===

積分区間には $f(x)=0$ の解を使い、$f(x)<0$ の区間では積分の符号を変える。

<code>
-3.63423216189401010467
-0.49319311896925860372
 0.41073330472197050467
 2.71669197614524948747
</code>

<file rsplus report4-2.R[enable_line_numbers="true"]>
# レポート
# シンプソンの公式
# 積分範囲を分ける方法

# 被積分関数の定義
f <- function(x) {
  return(-x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2)
}

simpson <- function(a, b) {
  # 幅
  h <- (b - a) / n
  
  s1 <- f(a + h * 0.5)
  s2 <- 0
  for (i in 1 : (n - 1)) {
    x <- a + i * h
    s1 <- s1 + f(x + h * 0.5)
    s2 <- s2 + f(x)
  }
  return((f(a) + 4 * s1 + 2 * s2 + f(b)) * h / 6)
}

# 区間数
n <- 1000

# f(x)=0 の解
x1 <- -3.63423216189401010467
x2 <- -0.49319311896925860372
x3 <- 0.41073330472197050467
x4 <- 2.71669197614524948747

s <- simpson(x1, x2) - simpson(x2, x3) + simpson(x3, x4)

print(sprintf("S = %f", s), quote = FALSE)
</file>

同じ処理を繰り返すために simpson() という関数を定義している。

----
==== 結果 ====

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

<code>
S = 87.251725
</code>

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

<code>
S = 138.341770
</code>

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

<code>
S = 246.731549
</code>

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

<code>
S = 447.533010
</code>

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

<code>
S = 795.972974
</code>