第2章: 行列式と連立一次方程式

A = np.array([[2, 1], [3, 4]]) B = np.array([[5, 6], [7, 8]]) det_A = np.linalg.det(A) det_B = np.linalg.det(B) det_AB = np.linalg.det(A @ B) print("行列式の性質の検証:") print(f"det(A) = {det_A:.4f}") print(f"det(B) = {det_B:.4f}") print(f"det(AB) = {det_AB:.4f}") print(f"det(A) × det(B) = {det_A * det_B:.4f}") print(f"det(AB) = det(A)×det(B)? {np.isclose(det_AB, det_A * det_B)}") # 転置 det_AT = np.linalg.det(A.T) print(f"\ndet(A^T) = {det_AT:.4f}") print(f"det(A) = det(A^T)? {np.isclose(det_A, det_AT)}")

# 連立方程式: 2x + 3y = 8, x - y = -1 A = np.array([[2, 3], [1, -1]]) b = np.array([8, -1]) # np.linalg.solve で解く x = np.linalg.solve(A, b) print("連立一次方程式の解:") print(f"2x + 3y = 8") print(f"x - y = -1") print(f"\n解: x = {x[0]:.4f}, y = {x[1]:.4f}") # 検算 b_check = A @ x print(f"\n検算: Ax = {b_check}") print(f"誤差: {np.linalg.norm(b - b_check):.2e}")

def cramers_rule(A, b): """ クラメルの公式による連立方程式の解法 x_i = det(A_i) / det(A) A_i: i列目をbで置き換えた行列 """ det_A = np.linalg.det(A) if np.abs(det_A) < 1e-10: raise ValueError("行列式がゼロのため解けません") n = len(b) x = np.zeros(n) for i in range(n): A_i = A.copy() A_i[:, i] = b # i列目をbで置き換え x[i] = np.linalg.det(A_i) / det_A return x # 同じ連立方程式を解く x_cramer = cramers_rule(A, b) print("クラメルの公式:") print(f"解: x = {x_cramer[0]:.4f}, y = {x_cramer[1]:.4f}") print(f"np.linalg.solveとの差: {np.linalg.norm(x - x_cramer):.2e}")

def gaussian_elimination(A, b): """ ガウスの消去法による連立方程式の解法 前進消去 → 後退代入 """ n = len(b) # 拡大係数行列を作成 Ab = np.hstack([A.astype(float), b.reshape(-1, 1)]) # 前進消去 for i in range(n): # ピボット選択（部分ピボット選択） max_row = np.argmax(np.abs(Ab[i:, i])) + i Ab[[i, max_row]] = Ab[[max_row, i]] # i行目を正規化 Ab[i] = Ab[i] / Ab[i, i] # i列目の下側をゼロにする for j in range(i+1, n): Ab[j] = Ab[j] - Ab[j, i] * Ab[i] # 後退代入 x = np.zeros(n) for i in range(n-1, -1, -1): x[i] = Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:]) return x # テスト A_test = np.array([[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]]) b_test = np.array([8.0, -11.0, -3.0]) x_gauss = gaussian_elimination(A_test.copy(), b_test.copy()) x_numpy = np.linalg.solve(A_test, b_test) print("ガウスの消去法:") print(f"ガウス法の解: {x_gauss}") print(f"NumPyの解: {x_numpy}") print(f"差: {np.linalg.norm(x_gauss - x_numpy):.2e}")

from scipy.linalg import lu # LU分解 A_lu = np.array([[4, 3], [6, 3]]) P, L, U = lu(A_lu) print("LU分解:") print(f"A =\n{A_lu}\n") print(f"P (置換行列) =\n{P}\n") print(f"L (下三角) =\n{L}\n") print(f"U (上三角) =\n{U}\n") # 検証: PA = LU print(f"PA =\n{P @ A_lu}\n") print(f"LU =\n{L @ U}\n") print(f"PA = LU? {np.allclose(P @ A_lu, L @ U)}") # LU分解を使った連立方程式の解法 b_lu = np.array([10, 12]) # Step 1: Ly = Pb を解く (前進代入) y = np.linalg.solve(L, P @ b_lu) # Step 2: Ux = y を解く (後退代入) x_lu = np.linalg.solve(U, y) print(f"\n連立方程式 Ax = b の解:") print(f"x = {x_lu}") # 検証 print(f"Ax = {A_lu @ x_lu}") print(f"b = {b_lu}")

# 階数の計算 A_rank = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) rank_A = np.linalg.matrix_rank(A_rank) print("階数の計算:") print(f"A =\n{A_rank}") print(f"rank(A) = {rank_A}") print(f"Aは3次正方行列だが rank < 3 なので特異行列\n") # 一意解を持つ例 A_unique = np.array([[1, 2], [3, 4]]) b_unique = np.array([5, 6]) Ab_unique = np.column_stack([A_unique, b_unique]) rank_A_u = np.linalg.matrix_rank(A_unique) rank_Ab_u = np.linalg.matrix_rank(Ab_unique) print("一意解の例:") print(f"rank(A) = {rank_A_u}") print(f"rank([A|b]) = {rank_Ab_u}") print(f"rank(A) = rank([A|b]) = 2 → 一意解") x_unique = np.linalg.solve(A_unique, b_unique) print(f"解: {x_unique}\n") # 解なしの例 A_no_sol = np.array([[1, 2], [2, 4]]) b_no_sol = np.array([3, 7]) # 第2式は第1式の2倍だが右辺が不整合 Ab_no_sol = np.column_stack([A_no_sol, b_no_sol]) rank_A_n = np.linalg.matrix_rank(A_no_sol) rank_Ab_n = np.linalg.matrix_rank(Ab_no_sol) print("解なしの例:") print(f"rank(A) = {rank_A_n}") print(f"rank([A|b]) = {rank_Ab_n}") print(f"rank(A) < rank([A|b]) → 解なし")

# 化学反応式: aFe + bO2 → cFe2O3 # 元素バランス: # Fe: a = 2c # O: 2b = 3c # これを行列形式で表現 # 係数行列（未知数を左辺、既知数を右辺に） # a - 2c = 0 # 2b - 3c = 0 # c = 1 (任意に固定) # 解きやすい形に変形 A_chem = np.array([[1, 0, -2], # Fe balance [0, 2, -3]]) # O balance # c=1とすると c = 1 b_chem = np.array([2*c, 3*c]) # 最小二乗解を求める x_chem = np.linalg.lstsq(A_chem[:, :2], b_chem, rcond=None)[0] a, b = x_chem print("化学反応式のバランシング:") print(f"Fe + O2 → Fe2O3") print(f"\n係数:") print(f"a (Fe) = {a:.1f}") print(f"b (O2) = {b:.1f}") print(f"c (Fe2O3) = {c:.1f}") print(f"\nバランスした反応式:") print(f"{int(a*2)}Fe + {int(b*2)}O2 → {int(c*2)}Fe2O3") # 検証 print(f"\n元素バランスの検証:") print(f"Fe: 左辺 = {int(a*2)}, 右辺 = {int(c*2)*2}") print(f"O: 左辺 = {int(b*2)*2}, 右辺 = {int(c*2)*3}")