Rで線形代数:逆行列

# --------------------------------------------------
# 逆行列のシミュレーション
# --------------------------------------------------

# --- シナリオ1: 基本的な逆行列の計算 ---

cat("--- シナリオ1: 基本的な逆行列の計算 ---\n")

# 逆行列が存在する行列 A を作成
# det(A) = 2*4 - 3*1 = 8 - 3 = 5 ≠ 0 なので逆行列は存在する
A <- matrix(c(2, 1, 3, 4), nrow = 2, ncol = 2)
cat("元の行列 A:\n")
print(A)

# Aの逆行列を計算
A_inv <- solve(A)
cat("\nAの逆行列 A_inv:\n")
print(A_inv)

# 確認: A %*% A_inv が単位行列 I になるか？
I <- A %*% A_inv
cat("\n確認 (A %*% A_inv):\n")
# 計算誤差を考慮して丸めて表示
print(round(I, 10))
cat("→ 単位行列になりました。\n\n")


# --- シナリオ2: 逆行列が存在しない場合（特異行列） ---

cat("--- シナリオ2: 逆行列が存在しない場合 ---\n")

# 逆行列が存在しない行列 B を作成
# 2行目が1行目の2倍になっており、線形独立ではない
# det(B) = 1*4 - 2*2 = 4 - 4 = 0
B <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2)
cat("特異行列 B:\n")
print(B)

# solve(B) を実行するとエラーになることを確認
# tryCatch でエラーを捕捉して、エラーメッセージを表示
result <- tryCatch(
  {
    solve(B)
  },
  error = function(e) {
    return(e) # エラーオブジェクトを返す
  }
)

cat("\n特異行列 B の逆行列を計算しようとすると...\n")
if (inherits(result, "error")) {
  cat("エラーが発生しました:\n")
  print(result$message)
  cat("→ 行列式が0のため、逆行列を計算できませんでした。\n\n")
}


# --- シナリオ3: 応用例（連立一次方程式の解法） ---

cat("--- シナリオ3: 応用例（連立一次方程式の解法） ---\n")

# 方程式:
# 2x + y = 8
#  x + 3y = 9

# 行列形式 Ax = b に変換
A_eq <- matrix(c(2, 1, 1, 3), nrow = 2, byrow = TRUE)
b_eq <- matrix(c(8, 9), nrow = 2)

cat("係数行列 A:\n")
print(A_eq)
cat("\n定数ベクトル b:\n")
print(b_eq)

# 解 x = A⁻¹b を計算
# 方法1: 逆行列を明示的に計算してから掛ける
A_eq_inv <- solve(A_eq)
x_sol1 <- A_eq_inv %*% b_eq

# 方法2: solve(A, b) を使う
x_sol2 <- solve(A_eq, b_eq)

cat("\n解 x (x_sol1):\n")
print(x_sol1)

cat("\n解 x (x_sol2):\n")
print(x_sol2)

# 検算: A %*% x が b に戻るか確認
check <- A_eq %*% x_sol2
cat("\n検算 (A %*% x):\n")
print(check)
cat("→ 元のベクトル b と一致しました。\n")

cat(paste0("\nよって、方程式の解は x = ", x_sol2[1], ", y = ", x_sol2[2], " です。\n"))