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

課題

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

求める面積

解答例

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

第2回課題の $f(x)=0$ の4つの解を使って積分範囲を3つに分割して積分を行う。 $f(x)<0$ の区間では積分の符号を変える。

report3.R
  2. # シンプソンの公式
  3.  
  4. # 被積分関数の定義
  5. f <- function(x) {
  6.   return(-x^4 - x^3 + 10 * x^2 + x - 2)
  7. }
  8.  
  9. # 区間数
  10. n <- 1000
  11.  
  12. # シンプソンの公式による [a,b] の数値積分
  13. simpson <- function(a, b) {
  14.   # 幅
  15.   h <- (b - a) / n
  16.  
  17.   s1 <- f(a + h * 0.5)
  18.   s2 <- 0
  19.   for (i in 1 : (n - 1)) {
  20.     x <- a + i * h
  21.     s1 <- s1 + f(x + h * 0.5)
  22.     s2 <- s2 + f(x)
  23.   }
  24.   s <- (f(a) + 4 * s1 + 2 * s2 + f(b)) * h / 6
  25.   return(s)
  26. }
  27.  
  28. x1 <- -3.63423216189423481381
  29. x2 <- -0.49319311896925860372
  30. x3 <- 0.41073330471801877684
  31. x4 <- 2.71669197614524948747
  32.  
  33. s1 <- simpson(x1, x2)
  34. s2 <- simpson(x2, x3)
  35. s3 <- simpson(x3, x4)
  36. print(sprintf("S1 = %f", s1), quote = FALSE)
  37. print(sprintf("S2 = %f", s2), quote = FALSE)
  38. print(sprintf("S3 = %f", s3), quote = FALSE)
  39. print(sprintf("S = %f", abs(s1) + abs(s2) + abs(s3)), quote = FALSE)

同じ処理を繰り返すために区間 $[a,b]$ の数値積分を行う simpson(a,b) という関数を定義している。 最後に負の面積を正にするため abs() を使用している。

結果

$f(x)=0$ の解やパラメーターの違いにより小数点以下の値が多少異なることがある。

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

S = 87.251731

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

S = 138.341753

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

S = 246.731532

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

S = 447.533001

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

S = 795.972953

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

S = 1372.366483