解线性方程组的直接方法：高斯消元法与其程序实现

1.顺序高斯消元法

设线性方程组$$\mathbf{A}\mathbf{x}=\mathbf{b}$$
如果$$a_{kk}^{\left(k\right)}\ne 0$$
可以通过高斯消元法转化为等价的三角形线性方程组：
$$\begin{gathered} \left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right] \\ \stackrel{\text{高斯消元法}}{\to }\left[\begin{array}{cccc}{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)}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & a_{nn}^{\left(n\right)}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1^{\left(1\right)}\\ b_2^{\left(2\right)}\\ ⋮\\ b_n^{\left(n\right)}\end{array}\right] \end{gathered}$$

2.列主元高斯消元法

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

• 当线性方程组的系数行列式$\left|\mathbf{A}\right|\ne 0$,且存在$1$个主元素为$0$，线性方程组有为一解，但是顺序高斯消元法不能解出。
• 小主元可以求解但是误差比较大。

通过交换行来避免小主元和$0$主元，然后继续使用消元法求解的方法，称为选主元高斯消元法
其改进之处在于：
选主元：
$$\left|a_{pk}^{(k)}\right|=\underset{k\le i\le n}{\max }\left|a_{ik}^{(k)}\right|$$

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

消元公式：
$$m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},i=k+1,\cdots ,n,$$
$$\{\begin{array}{ll}{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{array}$$

3.全主元高斯消元法

列主元高斯消元法虽然在一定程度上避免了零主元和小主元带来的问题，但它仅在当前列中选取主元。而全主元高斯消元法更为严格，它在整个未处理的子矩阵中选取主元，从而能进一步减少计算误差，提高求解的稳定性。

对于线性方程组 $\mathbf{A}\mathbf{x}=\mathbf{b}$，其中 $\mathbf{A}$ 是 $n\times n$ 系数矩阵，$\mathbf{b}$ 是 $n$ 维常数向量，求解步骤如下：
对于线性方程组 $\mathbf{A}\mathbf{x}=\mathbf{b}$，其中 $\mathbf{A}$ 是 $n\times n$ 系数矩阵，$\mathbf{b}$ 是 $n$ 维常数向量，求解步骤如下：

步骤 1：初始化

设初始矩阵为 ${\mathbf{A}}^{(1)}=\mathbf{A}$，${\mathbf{b}}^{(1)}=\mathbf{b}$，$k=1$。

步骤 2：选主元（在第 $k$ 步）

在子矩阵 ${\mathbf{A}}^{(k)}$ 的右下角 $(n-k+1)\times (n-k+1)$ 子矩阵中选取绝对值最大的元素作为主元，即
$$\left|a_{pq}^{(k)}\right|=\underset{k\le i\le n,k\le j\le n}{\max }\left|a_{ij}^{(k)}\right|$$
其中，$p$ 表示主元所在的行，$q$ 表示主元所在的列。

步骤 3：交换行和列

若 $p\ne k$，交换 ${\mathbf{A}}^{(k)}$ 的第 $k$ 行和第 $p$ 行，同时交换 ${\mathbf{b}}^{(k)}$ 的第 $k$ 个和第 $p$ 个元素；若 $q\ne k$，交换 ${\mathbf{A}}^{(k)}$ 的第 $k$ 列和第 $q$ 列。记录列交换的信息，用于后续恢复解的顺序。交换规则如下：
$$\{\begin{array}{ll}{a}_{kj}^{(k)}\leftrightarrow a_{pj}^{(k)},& j=1,\cdots ,n\\ b_k^{(k)}\leftrightarrow b_p^{(k)}& \end{array}$$
$$\{\begin{array}{ll}{a}_{ik}^{(k)}\leftrightarrow a_{iq}^{(k)},& i=1,\cdots ,n\end{array}$$

步骤 4：消元计算

计算乘数：
$$m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},i=k+1,\cdots ,n$$
更新系数矩阵和常数向量：
$$\{\begin{array}{ll}{a}_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ik}{a}_{kj}^{(k)},& i=k+1,\cdots ,n;j=k+1,\cdots ,n\\ b_i^{(k+1)}=b_i^{(k)}-m_{ik}{b}_k^{(k)},& i=k+1,\cdots ,n\end{array}$$

步骤 5：循环更新

令 $k=k+1$，若 $k<n$，返回步骤 2 继续进行；若 $k=n$，则消元过程结束。

步骤 6：回代求解

经过消元后，得到上三角方程组 ${\mathbf{A}}^{(n)}\mathbf{x}={\mathbf{b}}^{(n)}$，通过回代法求解 $\mathbf{x}$。回代公式为：
$$x_n=\frac{b_n^{(n)}}{a_{nn}^{(n)}}$$
$$x_i=\frac{b_i^{(i)}-{\sum }_{j=i+1}^n{a}_{ij}^{(i)}{x}_j}{a_{ii}^{(i)}},i=n-1,\cdots ,1$$

4.高斯——若尔当消元法

以下是高斯-若尔当消元法的分步解释和公式表示：

---

高斯-若尔当消元法（Gauss-Jordan Elimination）

核心思想：
在消元过程中，不仅通过前向消元将系数矩阵化为上三角矩阵（如普通高斯消元法），进一步通过反向消元将矩阵化为对角矩阵或单位矩阵，从而直接得到解，无需回代。

---

分步矩阵变换过程

初始线性方程组（增广矩阵形式）：
$$\left[\begin{array}{ccccc}{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)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}^{(1)}& a_{n2}^{(1)}& \cdots & a_{nn}^{(1)}& b_n^{(1)}\end{array}\right]$$

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

1. 选取主元：在第 $k$列中选取绝对值最大的元素作为主元（部分选主元法）。
2. 归一化主元行：将主元所在行除以主元值，使主元位置为 $1$。
3. 消去主元列的所有非主元元素：
   • 对 $i\ne k$ 的所有行，用主元行消去第 $k$ 列的元素，使其为$0$。

分步过程：
$$\begin{array}{rl}& \stackrel{\text{第1步消元}}{\to }\left[\begin{array}{ccccc}1& a_{12}^{(2)}& \cdots & a_{1n}^{(2)}& b_1^{(2)}\\ 0& a_{22}^{(2)}& \cdots & a_{2n}^{(2)}& b_2^{(2)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& a_{n2}^{(2)}& \cdots & a_{nn}^{(2)}& b_n^{(2)}\end{array}\right]\\ & \stackrel{\text{第2步消元}}{\to }\left[\begin{array}{ccccc}1& 0& \cdots & a_{1n}^{(3)}& b_1^{(3)}\\ 0& 1& \cdots & a_{2n}^{(3)}& b_2^{(3)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & a_{nn}^{(3)}& b_n^{(3)}\end{array}\right]\\ & ⋮\\ & \stackrel{\text{第n步消元}}{\to }\left[\begin{array}{ccccc}1& 0& \cdots & 0& b_1^{(n+1)}\\ 0& 1& \cdots & 0& b_2^{(n+1)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& b_n^{(n+1)}\end{array}\right]\end{array}$$
最终增广矩阵的最后一列即为解：
$$x_1=b_1^{(n+1)},x_2=b_2^{(n+1)},\dots ,x_n=b_n^{(n+1)}.$$

---

元素更新公式

在第 (k) 步消元中，假设主元已归一化为 1（即 $a_{kk}^{(k)}=1$），则对每行 $i\ne k$：
$$\{\begin{array}{ll}{a}_{ij}^{(k+1)}=a_{ij}^{(k)}-a_{ik}^{(k)}\cdot a_{kj}^{(k)},& \forall j\ge k,\\ b_i^{(k+1)}=b_i^{(k)}-a_{ik}^{(k)}\cdot b_k^{(k)},& \forall i\ne k.\end{array}$$

---

特点与应用

1. 直接求解：最终矩阵为单位矩阵，解直接显式给出，无需回代。
2. 计算复杂度：计算量约为 (O(n^3))，略高于普通高斯消元法。
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)