Rで逆2次補間法

# 必要なライブラリを読み込みます
library(ggplot2)
library(dplyr)

# -----------------------------------------------------------------
# 1. シミュレーションの対象となる関数を定義 (これまでと同じ)
# -----------------------------------------------------------------
# f(x) = x^3 - x - 1
f <- function(x) {
  return(x^3 - x - 1)
}

# -----------------------------------------------------------------
# 2. 逆2次補間法のシミュレーションを行う関数を定義
# -----------------------------------------------------------------
inverse_quadratic_simulation <- function(func, x0, x1, x2, tol = 1e-12, max_iter = 10) {
  cat("--- 逆2次補間法シミュレーション開始 ---\n\n")
  cat(sprintf("対象関数: x^3 - x - 1\n"))
  cat(sprintf("初期点: x0=%.1f, x1=%.1f, x2=%.1f\n", x0, x1, x2))
  cat(sprintf("許容誤差: %e\n\n", tol))

  # 各ステップの結果を保存するデータフレーム
  history <- data.frame(
    iter = integer(),
    x_new = double(),
    fx_new = double()
  )

  # 反復計算
  for (i in 1:max_iter) {
    # 現在の3点での関数値を計算
    y0 <- func(x0)
    y1 <- func(x1)
    y2 <- func(x2)

    # yの値が重複すると計算できないためチェック (簡易的)
    if (y0 == y1 || y0 == y2 || y1 == y2) {
      stop("関数値yが重複しました。異なる関数値を持つ3点を選択してください。")
    }

    # ラグランジュ補間による次の近似解の計算
    x_new <- (y1 * y2 / ((y0 - y1) * (y0 - y2))) * x0 +
      (y0 * y2 / ((y1 - y0) * (y1 - y2))) * x1 +
      (y0 * y1 / ((y2 - y0) * (y2 - y1))) * x2

    fx_new <- func(x_new)

    cat(sprintf("【反復 %d回目】\n", i))
    cat(sprintf("  次の近似解 x_new = %.12f\n", x_new))
    cat(sprintf("  その関数値 f(x_new) = %.2e\n", fx_new))

    history <- rbind(history, data.frame(
      iter = i,
      x_new = x_new,
      fx_new = fx_new
    ))

    # 収束判定
    if (abs(fx_new) < tol || abs(x_new - x2) < tol) {
      cat("\n--- 収束しました ---\n\n")
      break
    }

    # 次の反復のために点を更新
    x0 <- x1
    x1 <- x2
    x2 <- x_new
  }

  if (i == max_iter) {
    cat("\n--- 最大反復回数に達しました ---\n")
  }

  cat(sprintf("最終的な解の近似値 = %.12f\n", x_new))

  return(history)
}

# -----------------------------------------------------------------
# 3. シミュレーションの実行
# -----------------------------------------------------------------
# 3つの初期点を設定
sim_history_iqi <- inverse_quadratic_simulation(f, x0 = 2.0, x1 = 1.5, x2 = 1.0)

# -----------------------------------------------------------------
# 4. ggplot2による結果の可視化
# -----------------------------------------------------------------
# 最終的な解の値を取得
final_solution_value <- tail(sim_history_iqi, 1)$x_new

# 収束プロセスのプロット
iqi_process_plot <- ggplot(sim_history_iqi, aes(x = iter, y = x_new)) +
  geom_hline(yintercept = final_solution_value, linetype = "dashed", color = "red", linewidth = 1) +
  geom_text(aes(x = 0.5, y = final_solution_value, label = sprintf("真の解\n(約 %.6f)", final_solution_value)),
    color = "red", hjust = 0, vjust = -0.2
  ) +
  geom_line(color = "gray50", linewidth = 1) +
  geom_point(color = "#7CAE00", size = 5, stroke = 1.5) +
  labs(
    title = "逆2次補間法の収束プロセス",
    subtitle = "近似解がどのように真の解に収束するか (単独使用時)",
    x = "反復回数 (iter)",
    y = "近似解 (x_new) の値"
  ) +
  scale_x_continuous(breaks = 1:nrow(sim_history_iqi)) +
  theme_minimal(base_family = "sans")

print(iqi_process_plot)