Rで線形代数:クラメールの公式

Rで 線形代数:クラメールの公式 を確認します。

1. クラメールの公式（Cramer’s Rule）とは

目的

クラメールの公式は、連立一次方程式の解を、行列式を使って表現する公式です。 特に、解の存在や一意性に関する理論的な議論や、解がパラメータにどう依存するかを分析する際に有用です。

対象となる方程式

クラメールの公式が使えるのは、以下の条件を満たす連立一次方程式です。

  1. 方程式の数と未知数の数が同じであること。 （係数行列が n × n の正方行列になる）
  2. 係数行列 A の行列式が0でないこと (det(A) ≠ 0)。 （つまり、方程式の解がただ一つだけ存在する場合）

公式

n個の未知数 x₁, x₂, ..., xₙ を持つ以下の連立一次方程式を考えます。

a₁₁x₁ + a₁₂x₂ + ... + a₁ₙxₙ = b₁
a₂₁x₁ + a₂₂x₂ + ... + a₂ₙxₙ = b₂
...
aₙ₁x₁ + aₙ₂x₂ + ... + aₙₙxₙ = bₙ

これは行列で Ax = b と書けます。

このとき、各未知数 xᵢ は以下の公式で求められます。

xᵢ = det(Aᵢ) / det(A)

ここで、

  • det(A): 元の係数行列 A の行列式。
  • det(Aᵢ): 行列 A の i 番目の列を、定数項のベクトル b で置き換えた行列 Aᵢ の行列式。

具体例 (2×2の場合)

以下の連立方程式を考えます。

a x + b y = e
c x + d y = f

係数行列 A と定数ベクトル b は、

A = [[a, b], [c, d]]
b = [[e], [f]]

クラメールの公式によると、解 x と y は、

  • xを求める:

    A₁ = [[e, b], [f, d]] （Aの1列目をbで置き換え）

    x = det(A₁) / det(A) = (ed - bf) / (ad - bc)

  • yを求める:

    A₂ = [[a, e], [c, f]] （Aの2列目をbで置き換え）

    y = det(A₂) / det(A) = (af - ec) / (ad - bc)

となります。これは連立方程式の解の公式そのものです。

2. Rによるシミュレーション

Rを使って、クラメールの公式が正しく解を与えることを確認します。3×3の連立一次方程式を例にとり、公式に従って手計算し、Rの solve() 関数（効率的な内部アルゴリズムを使用）の結果と比較します。

Rコード

# --------------------------------------------------
# クラメールの公式のシミュレーション
# --------------------------------------------------

# --- 問題設定 ---
# 以下の3x3の連立一次方程式を解く
# 2x +  y -  z =  8
# -3x - y + 2z = -11
# -2x +  y + 2z = -3

# 行列形式 Ax = b に変換
A <- matrix(c(
  2, -3, -2,
  1, -1, 1,
  -1, 2, 2
), nrow = 3, byrow = FALSE)

b <- c(8, -11, -3)

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

# --- クラメールの公式で解を計算 ---

# Step 1: det(A) を計算
det_A <- det(A)

# det(A)が0でないことを確認 (0なら公式は使えない)
if (abs(det_A) < 1e-9) {
  stop("det(A)が0のため、クラメールの公式は適用できません。")
} else {
  cat("\ndet(A) =", det_A, "(≠ 0 なので、解は一意に存在します)\n")
}

# Step 2: det(A1), det(A2), det(A3) を計算
# A1: Aの1列目をbで置き換え
A1 <- A
A1[, 1] <- b
det_A1 <- det(A1)
cat("det(A1) =", det_A1, "\n")

# A2: Aの2列目をbで置き換え
A2 <- A
A2[, 2] <- b
det_A2 <- det(A2)
cat("det(A2) =", det_A2, "\n")

# A3: Aの3列目をbで置き換え
A3 <- A
A3[, 3] <- b
det_A3 <- det(A3)
cat("det(A3) =", det_A3, "\n")

# Step 3: 解を計算 x_i = det(A_i) / det(A)
x1 <- det_A1 / det_A
x2 <- det_A2 / det_A
x3 <- det_A3 / det_A

cramer_solution <- c(x1, x2, x3)
cat("\nクラメールの公式による解 (x1, x2, x3):\n")
print(cramer_solution)


# --- 比較: Rのsolve()関数による解法 ---

# solve(A, b) の内部でクラメールの公式は利用されていません
solve_solution <- solve(A, b)

cat("\nRのsolve()関数による解:\n")
print(solve_solution)

# --- 結果の検証 ---
cat("\n両者の解は一致するか？ ->", isTRUE(all.equal(cramer_solution, solve_solution)), "\n")

# 検算: A %*% x が b に戻るか確認
cat("\n検算 (A %*% x):\n")
check <- A %*% cramer_solution
print(check)
cat("→ 元のベクトル b と一致しました。\n")
係数行列 A:
     [,1] [,2] [,3]
[1,]    2    1   -1
[2,]   -3   -1    2
[3,]   -2    1    2

定数ベクトル b:
[1]   8 -11  -3

det(A) = -1 (≠ 0 なので、解は一意に存在します)
det(A1) = -2 
det(A2) = -3 
det(A3) = 1 

クラメールの公式による解 (x1, x2, x3):
[1]  2  3 -1

Rのsolve()関数による解:
[1]  2  3 -1

両者の解は一致するか？ -> TRUE 

検算 (A %*% x):
     [,1]
[1,]    8
[2,]  -11
[3,]   -3
→ 元のベクトル b と一致しました。

以上です。