====== 【数値情報処理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行目を変えること）

<file rsplus report5.R[enable_line_numbers="true",highlight_lines_extra="5,8,9,10,11,13,14,15,16,24,25,34"]>
# Romberg1段公式を使ったニュートン法

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

# 中心差分
f2 <- function(x, h) {
  return((f(x + h) - f(x - h)) / (2 * h))
}

# Romberg1段公式
fr <- function(x, h) {
  return((4 * f2(x, h) - f2(x, 2 * h)) / 3)
}

# 反復の回数
iter <- 100

# δ
delta <- 1e-10

# h
h <- 1e-1

# 初期値
x0 <- 0.5

f0 <- f(x0)

# 反復
for (k in 1 : iter) {
  x0 <- x0 - f0 / fr(x0, h)
  print(sprintf("%d回目 %.20f", k, x0), quote = FALSE)
  f0 <- f(x0)
  # |f(x)|<δ を満たせば終了
  if (abs(f0) < delta) {
    break
  }
}
</file>

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