Rでハレー法

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

# -----------------------------------------------------------------
# 1. 関数、1階導関数、2階導関数を定義
# -----------------------------------------------------------------
# f(x) = x^3 - x - 1
f <- function(x) {
  return(x^3 - x - 1)
}

# f'(x) = 3x^2 - 1
f_prime <- function(x) {
  return(3 * x^2 - 1)
}

# f''(x) = 6x
f_double_prime <- function(x) {
  return(6 * x)
}

# -----------------------------------------------------------------
# 2. ハレー法のシミュレーションを行う関数を定義
# -----------------------------------------------------------------
halley_simulation <- function(func, func_prime, func_double_prime, x0, tol = 1e-15, max_iter = 10) {
  cat("--- ハレー法シミュレーション開始 ---\n\n")
  cat(sprintf("対象関数: x^3 - x - 1\n"))
  cat(sprintf("初期値: x0 = %.1f\n", x0))
  cat(sprintf("許容誤差: %e\n\n", tol))

  x <- x0

  history <- data.frame(
    iter = integer(),
    x = double(),
    fx = double(),
    x_new = double(),
    fx_new = double()
  )

  for (i in 1:max_iter) {
    fx <- func(x)
    f_prime_x <- func_prime(x)
    f_double_prime_x <- func_double_prime(x)

    # 分母がゼロに近い場合は停止
    denominator <- 2 * f_prime_x^2 - fx * f_double_prime_x
    if (abs(denominator) < 1e-12) {
      cat("計算式の分母がゼロに近いため、計算を停止します。\n")
      break
    }

    # ハレー法の反復式
    x_new <- x - (2 * fx * f_prime_x) / denominator
    fx_new <- func(x_new)

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

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

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

  if (i == max_iter) {
    cat("\n--- 最大反復回数に達しました ---\n")
  }
  cat(sprintf("最終的な解の近似値 = %.15f\n", x))

  return(history)
}

# -----------------------------------------------------------------
# 3. シミュレーションの実行
# -----------------------------------------------------------------
sim_history_halley <- halley_simulation(f, f_prime, f_double_prime, x0 = 2.0)

# -----------------------------------------------------------------
# 4. ggplot2による結果の可視化 (2つのグラフを生成)
# -----------------------------------------------------------------
# --- グラフ1: 収束プロット ---
final_solution_value <- tail(sim_history_halley, 1)$x_new
halley_process_plot <- ggplot(sim_history_halley, 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 = "darkgreen", size = 5, stroke = 1.5) +
  labs(
    title = "ハレー法の収束プロセス",
    x = "反復回数 (iter)",
    y = "近似解 (x_new) の値"
  ) +
  scale_x_continuous(breaks = 1:nrow(sim_history_halley)) +
  theme_minimal(base_family = "sans")

# --- グラフ2: 求解プロセスの可視化 ---
curve_data <- data.frame(x = seq(1.2, 2.1, length.out = 200)) %>%
  mutate(y = f(x))
final_solution_point <- tail(sim_history_halley, 1)

visualization_plot <- ggplot() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_line(data = curve_data, aes(x = x, y = y), color = "blue", linewidth = 1.2) +

  # 各ステップの(x, f(x))から(x_new, 0)への移動を表現
  geom_segment(
    data = sim_history_halley,
    aes(x = x, y = fx, xend = x_new, yend = 0, color = iter),
    arrow = arrow(length = unit(0.2, "cm")), # 矢印で表現
    linewidth = 1
  ) +
  geom_point(data = sim_history_halley, aes(x = x, y = fx, color = iter), size = 4) +
  geom_point(aes(x = final_solution_point$x_new, y = 0), size = 5, shape = 8, color = "red") +
  geom_text(
    aes(
      x = final_solution_point$x_new, y = 0.5,
      label = sprintf("近似解\nx = %.5f", final_solution_point$x_new)
    ),
    color = "red", hjust = 0, vjust = 0
  ) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen") +
  labs(
    title = "ハレー法による求解プロセス",
    x = "x", y = "f(x)"
  ) +
  theme_minimal(base_family = "sans") +
  theme(legend.position = "none")

# 2つのプロットを表示
print(halley_process_plot)
print(visualization_plot)