Rで二分法

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

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

# -----------------------------------------------------------------
# 2. 二分法のシミュレーションを行う関数を定義
# -----------------------------------------------------------------
bisection_simulation <- function(func, a, b, tol = 1e-6, max_iter = 20) {
  # f(a)とf(b)の符号が同じ場合はエラー
  if (func(a) * func(b) >= 0) {
    stop("初期区間[a, b]でf(a)とf(b)の符号が同じです。異なる符号になるように設定してください。")
  }

  cat("--- 二分法シミュレーション開始 ---\n\n")
  cat(sprintf("対象関数: x^3 - x - 1\n"))
  cat(sprintf("初期区間: [%.1f, %.1f]\n", a, b))
  cat(sprintf("許容誤差: %e\n\n", tol))

  # 各ステップの結果を保存するためのデータフレームを初期化
  history <- data.frame(
    iter = integer(),
    a = double(),
    b = double(),
    c = double(),
    fa = double(),
    fb = double(),
    fc = double(),
    width = double()
  )

  # 反復計算
  for (i in 1:max_iter) {
    c <- (a + b) / 2 # 中点を計算
    fa <- func(a)
    fb <- func(b)
    fc <- func(c)

    # 現在のステップの情報をデータフレームに追加
    history <- rbind(history, data.frame(
      iter = i, a = a, b = b, c = c,
      fa = fa, fb = fb, fc = fc, width = b - a
    ))

    # 現在の状態を出力
    cat(sprintf("【反復 %d回目】\n", i))
    cat(sprintf("  区間 = [%.6f, %.6f], 区間の幅 = %.6f\n", a, b, b - a))
    cat(sprintf("  中点 c = %.6f\n", c))
    cat(sprintf("  f(a) = %.6f, f(c) = %.6f\n", fa, fc))

    # 収束判定
    if ((b - a) / 2 < tol || fc == 0) {
      cat("\n--- 収束しました ---\n\n")
      break
    }

    # 区間の更新
    if (fa * fc < 0) {
      b <- c # 解は [a, c] にある
      cat("  f(a)とf(c)の符号が異なるため、次の区間は [a, c] になります。\n\n")
    } else {
      a <- c # 解は [c, b] にある
      cat("  f(a)とf(c)の符号が同じため、次の区間は [c, b] になります。\n\n")
    }
  }

  if (i == max_iter && (b - a) / 2 >= tol) {
    cat("\n--- 最大反復回数に達しました ---\n")
  }

  cat(sprintf("最終的な解の近似値 c = %.6f\n", c))
  cat(sprintf("そのときの関数値 f(c) = %e\n", fc))

  return(history)
}

# -----------------------------------------------------------------
# 3. シミュレーションの実行
# -----------------------------------------------------------------
# 初期値を設定
a_init <- 1.0
b_init <- 2.0
tolerance <- 1e-5

# シミュレーションを実行し、計算過程のデータを取得
simulation_history <- bisection_simulation(f, a_init, b_init, tol = tolerance)


# -----------------------------------------------------------------
# 4. ggplot2による結果の可視化
# -----------------------------------------------------------------
# 関数の曲線を描画するためのデータを作成
curve_data <- data.frame(x = seq(0.8, 2.2, length.out = 200)) %>%
  mutate(y = f(x))

# 最終的な解
final_solution <- tail(simulation_history, 1)

# プロットを作成
bisection_plot <- ggplot() +
  # y=0 の線 (解の目標)
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +

  # f(x)のグラフ
  geom_line(data = curve_data, aes(x = x, y = y), color = "blue", linewidth = 0.1) +

  # 各反復ステップでの中点 (c, f(c)) をプロット
  geom_point(data = simulation_history, aes(x = c, y = fc, color = iter), size = 3.0) +

  # 各反復ステップでの探索区間 [a, b] をy軸の下部に表示
  geom_segment(
    data = simulation_history,
    aes(x = a, xend = b, y = -1.5 - iter * 0.4, yend = -1.5 - iter * 0.4, color = iter),
    linewidth = 1.5
  ) +

  # テキストをセグメントの右横に表示
  geom_text(
    data = simulation_history,
    aes(x = b + 0.02, y = -1.5 - iter * 0.4, label = paste0(iter, " 回目")),
    size = 3, hjust = 0, vjust = 0.5, color = "black"
  ) +

  # 最終的な解を点とテキストで強調
  geom_point(data = final_solution, aes(x = c, y = 0), size = 5, shape = 8, color = "red") +
  geom_text(
    data = final_solution,
    aes(x = c, y = 0.5, label = sprintf("近似解\nx = %.5f", c)),
    color = "red", hjust = 0, vjust = 0
  ) +

  # スケールとラベルの設定
  scale_color_gradient(low = "lightblue", high = "darkred") +
  coord_cartesian(ylim = c(-9, 6)) +
  labs(
    title = "二分法による方程式 f(x) = x^3 - x - 1 の求解プロセス",
    subtitle = "区間が徐々に狭まり、解に収束していく様子",
    x = "x",
    y = "f(x)"
  ) +
  theme_minimal(base_family = "sans") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )


# プロットを表示
print(bisection_plot)