Rで理論化学:電離平衡

# --- (2) 濃度を変化させたときの電離度の変化 ---

# C0 を 10^-1 ~ 10^-6 mol/L まで対数的に100点とる(グラフ描画用)
C0_seq <- 10^seq(-1, -6, length.out = 100)

# 各濃度について厳密解の alpha を計算する関数
calc_alpha_exact <- function(C0, Ka) {
  (-Ka + sqrt(Ka^2 + 4 * C0 * Ka)) / (2 * C0)
}

alpha_seq <- calc_alpha_exact(C0_seq, Ka)

# 近似解(比較用)。C0が小さくなるとalphaが1に近づき近似が破綻することを見る
alpha_seq_approx <- sqrt(Ka / C0_seq)
alpha_seq_approx[alpha_seq_approx > 1] <- NA # 物理的に意味を持たない領域は除外

# データフレームにまとめる
df <- data.frame(
  C0 = C0_seq,
  alpha_exact = alpha_seq,
  alpha_approx = alpha_seq_approx
)

# ggplot2で片対数プロット(x軸: 濃度の対数、y軸: 電離度)
library(ggplot2)

df_long <- data.frame(
  C0 = rep(df$C0, 2),
  alpha = c(df$alpha_exact, df$alpha_approx),
  type = rep(c("厳密解", "近似解"), each = nrow(df))
)

p <- ggplot(df_long, aes(x = C0, y = alpha, color = type, linetype = type)) +
  geom_line(linewidth = 1) +
  scale_x_log10(
    labels = scales::trans_format("log10", scales::math_format(10^.x))
  ) +
  labs(
    title = "酢酸の電離度と濃度の関係(オストワルトの希釈律)",
    x = expression(C[0] ~ "(mol/L, 対数軸)"),
    y = expression(alpha ~ "(電離度)"),
    color = NULL, linetype = NULL
  ) +
  theme_minimal(base_size = 13)

print(p)

# オストワルトの希釈律の定量的確認:
C0_check <- c(1e-2, 1e-4)
alpha_check_exact <- calc_alpha_exact(C0_check, Ka)
alpha_check_approx <- sqrt(Ka / C0_check)

ratio_C <- C0_check[1] / C0_check[2]
ratio_alpha_exact <- alpha_check_exact[2] / alpha_check_exact[1]
ratio_alpha_approx <- alpha_check_approx[2] / alpha_check_approx[1]

cat("=== オストワルトの希釈律の確認 ===\n")
cat(sprintf("濃度比 C0(1e-2)/C0(1e-4) = %.1f\n", ratio_C))
cat(sprintf("電離度比(厳密解) alpha(1e-4)/alpha(1e-2) = %.2f\n", ratio_alpha_exact))
cat(sprintf(
  "電離度比(近似解, sqrt則) = %.2f (定義上 sqrt(%.0f) = %.2f と厳密に一致)\n",
  ratio_alpha_approx, ratio_C, sqrt(ratio_C)
))
cat(sprintf(
  "両者の差 = %.2f  ← (1-alpha)を1とみなす近似がどれだけ破れているかの目安\n",
  ratio_alpha_approx - ratio_alpha_exact
))

# 参考: それぞれの濃度でのalpha自体も表示しておく
cat("\n--- 各濃度での電離度(厳密解) ---\n")
cat(sprintf("C0 = %.0e : alpha = %.4f\n", C0_check[1], alpha_check_exact[1]))
cat(sprintf("C0 = %.0e : alpha = %.4f\n", C0_check[2], alpha_check_exact[2]))