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

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

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

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

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

積分区間には $f(x)=0$ の解のうち、端の二つを使う。

-3.63423216189401010467
 2.71669197614524948747

report4-1.R

# レポート
# シンプソンの公式
# 絶対値を使う方法
 
# 被積分関数の定義
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)

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

-3.63423216189401010467
-0.49319311896925860372
 0.41073330472197050467
 2.71669197614524948747

report4-2.R

# レポート
# シンプソンの公式
# 積分範囲を分ける方法
 
# 被積分関数の定義
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)

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

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

S = 87.251725

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

S = 138.341770

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

S = 246.731549

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

S = 447.533010

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

S = 795.972974