Rではさみうち法

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

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

# -----------------------------------------------------------------
# 2. はさみうち法のシミュレーションを行う関数を定義
# -----------------------------------------------------------------
false_position_simulation <- function(func, a, b, tol = 1e-6, max_iter = 20) {
  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()
  )

  # 反復計算
  for (i in 1:max_iter) {
    fa <- func(a)
    fb <- func(b)

    # 次の近似解cを計算 (はさみうち法の核)
    c <- (a * fb - b * fa) / (fb - fa)
    fc <- func(c)

    history <- rbind(history, data.frame(
      iter = i, a = a, b = b, c = c,
      fa = fa, fb = fb, fc = fc
    ))

    cat(sprintf("【反復 %d回目】\n", i))
    cat(sprintf("  区間 = [%.6f, %.6f]\n", a, b))
    cat(sprintf("  次の点 c = %.6f\n", c))
    cat(sprintf("  f(a) = %.6f, f(c) = %.6f\n", fa, fc))

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

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

  if (i == max_iter && abs(fc) >= 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_fp <- false_position_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_fp <- tail(simulation_history_fp, 1)

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

  # 各反復ステップでの弦 (a,f(a))-(b,f(b)) を描画
  geom_segment(
    data = simulation_history_fp,
    aes(x = a, y = fa, xend = b, yend = fb, color = iter),
    linetype = "dotted", linewidth = 0.8
  ) +

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

  # 各反復ステップで計算された点 c をx軸上にプロット
  geom_point(data = simulation_history_fp, aes(x = c, y = 0, color = iter), size = 3) +

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

  # スケールとラベルの設定
  scale_color_gradient(low = "lightblue", high = "darkred") +
  labs(
    title = "はさみうち法による方程式 f(x) = x^3 - x - 1 の求解プロセス",
    subtitle = "弦とx軸の交点を次の探索点とし、解に収束していく様子",
    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(false_position_plot)