Rでニュートン法

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

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

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


# -----------------------------------------------------------------
# 2. ニュートン法のシミュレーションを行う関数を定義
# -----------------------------------------------------------------
newton_simulation <- function(func, func_prime, x0, tol = 1e-12, 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

  # 開始点(x)も保存するように列を追加
  history <- data.frame(
    iter = integer(),
    x = double(), # 開始点
    fx = double(), # f(開始点)
    x_new = double(), # 次の点
    fx_new = double()
  )

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

    if (abs(f_prime_x) < 1e-12) {
      cat("導関数f'(x)の値がゼロに近いため、計算を停止します。\n")
      break
    }

    x_new <- x - fx / f_prime_x
    fx_new <- func(x_new)

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

    # 開始点(x, fx)もデータフレームに保存
    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("最終的な解の近似値 = %.12f\n", x))

  return(history)
}

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

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

# プロットを作成
tangent_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 = 1.2) +

  # 各反復ステップでの接線を描画
  # (x, f(x)) から (x_new, 0) へ線を引く
  geom_segment(
    data = sim_history_newton,
    aes(x = x, y = fx, xend = x_new, yend = 0, color = iter),
    linetype = "dashed", linewidth = 1
  ) +

  # 接点をプロット
  geom_point(data = sim_history_newton, aes(x = x, y = fx, color = iter), size = 4) +
  # x軸との交点をプロット
  geom_point(data = sim_history_newton, aes(x = x_new, y = 0, color = iter), size = 4, shape = 1) +

  # 最終的な解を強調
  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 = "lightblue", high = "purple4") +
  labs(
    title = "ニュートン法による求解プロセス（接線による可視化）",
    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" # 凡例は非表示
  )

# 2つのプロットを表示
print(newton_process_plot)
print(tangent_plot)