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

課題

$f(x)=0$ の解をニュートン法で求めるRのプログラムを第2回のプログラムを元に作って提出しなさい。

ただし、ニュートン法で使う $f(x)$ の微分は $h=0.1$ とした Romberg 1段公式による数値微分を使うこと。

解答例

$f(x)=-x^{4}-x^{3}+10x^{2}+x-2$ の場合（他の場合は 5行目を変えること）

report5.R
  2. # Romberg1段公式を使ったニュートン法
  3.  
  4. # 関数の定義
  5. f <- function(x) {
  6.   return(-x ^ 4 - x ^ 3 + 10 * x ^ 2 + x - 2)
  7. }
  8.  
  9. # 中心差分
  10. f2 <- function(x, h) {
  11.   return((f(x + h) - f(x - h)) / (2 * h))
  12. }
  13.  
  14. # Romberg1段公式
  15. fr <- function(x, h) {
  16.   return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
  17. }
  18.  
  19. # 反復の回数
  20. iter <- 100
  21.  
  22. # δ
  23. delta <- 1e-10
  24.  
  25. # h
  26. h <- 1e-1
  27.  
  28. # 初期値
  29. x0 <- 0.5
  30.  
  31. f0 <- f(x0)
  32.  
  33. # 反復
  34. for (k in 1 : iter) {
  35.   x0 <- x0 - f0 / fr(x0, h)
  36.   print(sprintf("%d回目 %.20f", k, x0), quote = FALSE)
  37.   f0 <- f(x0)
  38.   # |f(x)|<δ を満たせば終了
  39.   if (abs(f0) < delta) {
  40.     break
  41.   }
  42. }

第2回のニュートン法で微分を fp という関数で定義していたが、そこを Romberg 1段公式に変えている（8〜16行目）。