解線性方程組的直接方法：高斯消元法與其程序實現

1.順序高斯消元法

設線性方程組 $\boldsymbol{Ax}=\boldsymbol{b}$
如果 $a_{kk}^{\left( k \right)}\ne 0$
可以通過高斯消元法轉化為等價的三角形線性方程組：
$\left[ \begin{matrix} a_{11}& a_{12}& \cdots& a_{1n}\\ a_{21}& a_{22}& \cdots& a_{2n}\\ \vdots& \vdots& \ddots& \vdots\\ a_{n1}& a_{n2}& \cdots& a_{nn}\\ \end{matrix} \right] \left[ \begin{array}{c} x_1\\ x_2\\ \vdots\\ x_n\\ \end{array} \right] =\left[ \begin{array}{c} b_1\\ b_2\\ \vdots\\ b_n\\ \end{array} \right] \\ \overset{\text{高斯消元法}}{\rightarrow}\left[ \begin{matrix} a_{11}^{\left( 1 \right)}& a_{12}^{\left( 1 \right)}& \cdots& a_{1n}^{\left( 1 \right)}\\ 0& a_{22}^{\left( 2 \right)}& \cdots& a_{2n}^{\left( 2 \right)}\\ \vdots& \vdots& \ddots& \vdots\\ 0& 0& \cdots& a_{nn}^{\left( n \right)}\\ \end{matrix} \right] \left[ \begin{array}{c} x_1\\ x_2\\ \vdots\\ x_n\\ \end{array} \right] =\left[ \begin{array}{c} b_{1}^{\left( 1 \right)}\\ b_{2}^{\left( 2 \right)}\\ \vdots\\ b_{n}^{\left( n \right)}\\ \end{array} \right]$

2.列主元高斯消元法

由于順序高斯消元法存在以下缺陷：

當線性方程組的系數行列式 $\left|\boldsymbol{A}\right| \ne 0$ ,且存在 $1$ 個主元素為 $0$ ，線性方程組有為一解，但是順序高斯消元法不能解出。
小主元可以求解但是誤差比較大。

通過交換行來避免小主元和 $0$ 主元，然后繼續使用消元法求解的方法，稱為選主元高斯消元法
其改進之處在于：
選主元：
$\left|a_{pk}^{(k)}\right| = \max_{k \le i \le n} \left|a_{ik}^{(k)}\right|$

交換行：
$\begin{cases} a_{kj}^{(k)} \leftrightarrow a_{pj}^{(k)}, & j = k, k+1, \cdots, n, \\ b_k^{(k)} \leftrightarrow b_p^{(k)}. \end{cases}$

消元公式：
$m_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}, \quad i = k+1, \cdots, n,$
$\begin{cases} a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik}a_{kj}^{(k)}, & i,j = k+1, \cdots, n, \\ b_i^{(k+1)} = b_i^{(k)} - m_{ik}b_k^{(k)}, & i = k+1, \cdots, n. \end{cases}$

3.全主元高斯消元法

列主元高斯消元法雖然在一定程度上避免了零主元和小主元帶來的問題，但它僅在當前列中選取主元。而全主元高斯消元法更為嚴格，它在整個未處理的子矩陣中選取主元，從而能進一步減少計算誤差，提高求解的穩定性。

對于線性方程組 $\boldsymbol{Ax}=\boldsymbol{b}$ ，其中 $\boldsymbol{A}$ 是 $n\times n$ 系數矩陣， $\boldsymbol{b}$ 是 $n$ 維常數向量，求解步驟如下：
對于線性方程組 $\boldsymbol{Ax}=\boldsymbol{b}$ ，其中 $\boldsymbol{A}$ 是 $n\times n$ 系數矩陣， $\boldsymbol{b}$ 是 $n$ 維常數向量，求解步驟如下：

步驟 1：初始化

設初始矩陣為 $\boldsymbol{A}^{(1)}=\boldsymbol{A}$ ， $\boldsymbol{b}^{(1)}=\boldsymbol{b}$ ， $k = 1$ 。

步驟 2：選主元（在第 $k$ 步）

在子矩陣 $\boldsymbol{A}^{(k)}$ 的右下角 $k+1)\times(n - k + 1)$ 子矩陣中選取絕對值最大的元素作為主元，即
$\left|a_{p q}^{(k)}\right|=\max_{k\leq i\leq n,k\leq j\leq n}\left|a_{i j}^{(k)}\right|$
其中， $p$ 表示主元所在的行， $q$ 表示主元所在的列。

步驟 3：交換行和列

若 $p\neq k$ ，交換 $\boldsymbol{A}^{(k)}$ 的第 $k$ 行和第 $p$ 行，同時交換 $\boldsymbol{b}^{(k)}$ 的第 $k$ 個和第 $p$ 個元素；若 $q\neq k$ ，交換 $\boldsymbol{A}^{(k)}$ 的第 $k$ 列和第 $q$ 列。記錄列交換的信息，用于后續恢復解的順序。交換規則如下：
$\begin{cases} a_{k j}^{(k)}\leftrightarrow a_{p j}^{(k)},&j = 1,\cdots,n\\ b_{k}^{(k)}\leftrightarrow b_{p}^{(k)} \end{cases}$
$\begin{cases} a_{i k}^{(k)}\leftrightarrow a_{i q}^{(k)},&i = 1,\cdots,n \end{cases}$

步驟 4：消元計算

計算乘數：
$m_{i k}=\frac{a_{i k}^{(k)}}{a_{k k}^{(k)}}, \quad i=k + 1,\cdots,n$
更新系數矩陣和常數向量：
$\begin{cases} a_{i j}^{(k + 1)}=a_{i j}^{(k)}-m_{i k}a_{k j}^{(k)},&i = k + 1,\cdots,n;j = k+1,\cdots,n\\ b_{i}^{(k + 1)}=b_{i}^{(k)}-m_{i k}b_{k}^{(k)},&i = k + 1,\cdots,n \end{cases}$

步驟 5：循環更新

令 $k = k + 1$ ，若 $k\lt n$ ，返回步驟 2 繼續進行；若 $k = n$ ，則消元過程結束。

步驟 6：回代求解

經過消元后，得到上三角方程組 $\boldsymbol{A}^{(n)}\boldsymbol{x}=\boldsymbol{b}^{(n)}$ ，通過回代法求解 $\boldsymbol{x}$ 。回代公式為：
$x_{n}=\frac{b_{n}^{(n)}}{a_{n n}^{(n)}}$
$x_{i}=\frac{b_{i}^{(i)}-\sum_{j = i+1}^{n}a_{i j}^{(i)}x_{j}}{a_{i i}^{(i)}}, \quad i=n - 1,\cdots,1$

4.高斯——若爾當消元法

以下是高斯-若爾當消元法的分步解釋和公式表示：

高斯-若爾當消元法（Gauss-Jordan Elimination）

核心思想：
在消元過程中，不僅通過前向消元將系數矩陣化為上三角矩陣（如普通高斯消元法），進一步通過反向消元將矩陣化為對角矩陣或單位矩陣，從而直接得到解，無需回代。

分步矩陣變換過程

初始線性方程組（增廣矩陣形式）：
$\left[\begin{array}{cccc|c} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)}\\ a_{21}^{(1)} & a_{22}^{(1)} & \cdots & a_{2n}^{(1)} & b_2^{(1)}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1}^{(1)} & a_{n2}^{(1)} & \cdots & a_{nn}^{(1)} & b_n^{(1)} \end{array}\right]$

第 $k$ 步消元（共 $n$ 步）：

選取主元：在第 $k$ 列中選取絕對值最大的元素作為主元（部分選主元法）。
歸一化主元行：將主元所在行除以主元值，使主元位置為 $1$ 。
消去主元列的所有非主元元素：
- 對 $\neq k$ 的所有行，用主元行消去第 $k$ 列的元素，使其為 $0$ 。

分步過程：
$\begin{aligned} &\xrightarrow{\text{第1步消元}} \begin{bmatrix} 1 & a_{12}^{(2)} & \cdots & a_{1n}^{(2)} & b_1^{(2)} \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} & b_2^{(2)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)} & b_n^{(2)} \end{bmatrix} \\[2ex] &\xrightarrow{\text{第2步消元}} \begin{bmatrix} 1 & 0 & \cdots & a_{1n}^{(3)} & b_1^{(3)} \\ 0 & 1 & \cdots & a_{2n}^{(3)} & b_2^{(3)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & a_{nn}^{(3)} & b_n^{(3)} \end{bmatrix} \\[2ex] &\quad\ \ \vdots \\[2ex] &\xrightarrow{\text{第n步消元}} \begin{bmatrix} 1 & 0 & \cdots & 0 & b_1^{(n+1)} \\ 0 & 1 & \cdots & 0 & b_2^{(n+1)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & b_n^{(n+1)} \end{bmatrix} \end{aligned}$
最終增廣矩陣的最后一列即為解：
$x_1 = b_1^{(n+1)}, \quad x_2 = b_2^{(n+1)}, \quad \dots, \quad x_n = b_n^{(n+1)}.$

元素更新公式

在第 (k) 步消元中，假設主元已歸一化為 1（即 $a_{kk}^{(k)} = 1$ ），則對每行 $\neq k$ ：
$\begin{cases} a_{ij}^{(k+1)} = a_{ij}^{(k)} - a_{ik}^{(k)} \cdot a_{kj}^{(k)}, & \forall j \geq k, \\ b_i^{(k+1)} = b_i^{(k)} - a_{ik}^{(k)} \cdot b_k^{(k)}, & \forall i \neq k. \end{cases}$

特點與應用

直接求解：最終矩陣為單位矩陣，解直接顯式給出，無需回代。
計算復雜度：計算量約為 (O(n^3))，略高于普通高斯消元法。
應用場景：
- 求矩陣的逆：對單位矩陣進行同步行變換。
- 簡化方程組：直接得到最簡形式（如線性代數中的行最簡形）。

一道實訓題

import numpy as npclass GaussianElimination:def __init__(self, A, b, sol_method="principal"):self.A = np.asarray(A, dtype=np.float64)self.b = np.asarray(b, dtype=np.float64)if len(self.b) != self.A.shape[0]:raise ValueError("右端向量維度與系數矩陣維度不匹配！")if self.A.shape[0] != self.A.shape[1]:raise ValueError("系數矩陣不是方陣，不能用高斯消元法求解！")self.n = self.A.shape[0]self.augmented_matrix = np.c_[self.A, self.b]self.x = Noneself.eps = Noneself.sol_method = sol_methodself.jordan_inverse_matrix = Noneif self.sol_method == "sequential":self._solve_sequential_()elif self.sol_method == "column":self._solve_column_pivot_element_()elif self.sol_method == "complete":self._solve_complete_pivot_element_()elif self.sol_method == "jordan":self._solve_jordan_()else:raise ValueError("僅支持(sequential, column, complete, jordan)")def _elimination_process(self, i, k):if self.augmented_matrix[k, k] == 0:raise ValueError("系數矩陣不滿足順序高斯消元法求解！")multiplier = self.augmented_matrix[i, k] / self.augmented_matrix[k, k]self.augmented_matrix[i, k:] -= multiplier * self.augmented_matrix[k, k:]def _back_substitution_process_(self):x = np.zeros(self.n)for k in range(self.n - 1, -1, -1):sum_ = np.dot(self.augmented_matrix[k, k + 1:self.n], x[k + 1:self.n])x[k] = (self.augmented_matrix[k, -1] - sum_) / self.augmented_matrix[k, k]return xdef _solve_sequential_(self):for k in range(self.n - 1):for i in range(k + 1, self.n):self._elimination_process(i, k)self.x = self._back_substitution_process_()self.eps = np.dot(self.A, self.x) - self.bdef _solve_column_pivot_element_(self):for k in range(self.n - 1):idx = np.argmax(np.abs(self.augmented_matrix[k:, k])) + k# 修正索引錯誤，直接用 idxif idx != k:self.augmented_matrix[[k, idx]] = self.augmented_matrix[[idx, k]]for i in range(k + 1, self.n):self._elimination_process(i, k)self.x = self._back_substitution_process_()self.eps = np.dot(self.A, self.x) - self.bdef _solve_complete_pivot_element_(self):self.x = np.zeros(self.n)column_index = np.linspace(0, self.n - 1, self.n, dtype=np.int64)for k in range(self.n - 1):max_x = np.max(np.abs(self.augmented_matrix[k:, k:k + 1]))id_r, id_c = np.where(np.abs(self.augmented_matrix[k:, k:k + 1]) == max_x)id_r, id_c = int(id_r[0]), int(id_c[0])id_r = id_r + kif id_r != k:self.augmented_matrix[[k, id_r]] = self.augmented_matrix[[id_r, k]]id_c = id_c + kif id_c != k:pos = column_index[id_c]column_index[id_c] = column_index[k]column_index[k] = posself.augmented_matrix[:, [k, id_c]] = self.augmented_matrix[:, [id_c, k]]for i in range(k + 1, self.n):self._elimination_process(i, k)solve_x = self._back_substitution_process_()for k in range(self.n):for j in range(self.n):if k == column_index[j]:self.x[k] = solve_x[j]breakself.eps = np.dot(self.A, self.x) - self.bdef _solve_jordan_(self):self.augmented_matrix = np.c_[self.augmented_matrix, np.eye(self.n)]for k in range(self.n):idx = np.argmax(np.abs(self.augmented_matrix[k:, k])) + kif idx != k:self.augmented_matrix[[k, idx]] = self.augmented_matrix[[idx, k]]if self.augmented_matrix[k, k] == 0:raise ValueError("系數矩陣不滿足高斯 - 若爾當消元法")self.augmented_matrix[k, :] /= self.augmented_matrix[k, k]for i in range(self.n):if i != k:factor = self.augmented_matrix[i, k]self.augmented_matrix[i, :] -= factor * self.augmented_matrix[k, :]self.jordan_inverse_matrix = self.augmented_matrix[:, self.n + 1:]self.x = self.augmented_matrix[:, -self.n - 1]self.eps = np.dot(self.A, self.x) - self.b# 給定的線性方程組
A = [[10, -7, 0, 1],[-3, 2.099999, 6, 2],[5, -1, 5, -1],[2, 1, 0, 2]
]
b = [8, 5.900001, 5, 1]# 使用列主元高斯消去法求解
solver = GaussianElimination(A, b, sol_method="column")# 輸出結果
print("方程組的解：", solver.x)
print("解的精度驗證：", solver.eps)