# 数値微分（前進差分,中央差分,Romberg1段公式）の誤差

# 関数
f <- function(x) {
  return(x ^ 5 - 4 * x ^ 4)
}

# 導関数（誤差確認用）
fp <- function(x) {
  return(5 * x ^ 4 - 16 * x ^ 3)
}

# 前進差分
f1 <- function(x, h) {
  return((f(x + h) - f(x)) / h)
}

# 中心差分
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)
}

x <- 0.1
h <- 10^(0 : -20)
e1 <- abs(fp(x) - f1(x, h))
e2 <- abs(fp(x) - f2(x, h))
er <- abs(fp(x) - fr(x, h))
range <- c(1e-16, 1)
plot(h, e1, ylab = "error", ylim = range, log = "xy", type = "o", col = "blue")
par(new = TRUE)
plot(h, e2, ylim = range, log = "xy", type = "o", col = "red",
     ann = FALSE, axes = FALSE)
par(new = TRUE)
plot(h, er, ylim = range, log = "xy", type = "o", col = "green",
     ann = FALSE, axes = FALSE)