高斯-约旦消元法

创建时间：2025-11-04

---

高斯-约旦消元法是线性代数中常用的工具，可以以 \(O(n^2m)\) 的时间复杂度求解线性方程组，其中 \(n\) 是未知数个数，\(m\) 是方程个数。

以P2455 [SDOI2006] 线性方程组为例。题目给定含 \(n\) 个方程的线性方程组：

\[ \begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \ \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}\]

判断该方程组是否有解。若有解是否有唯一解，判断是否有唯一解；若有唯一解，输出这组解。

注意到这个方程组实际上只与 \(a\) 和 \(b\) 有关，故可以以矩阵的形式更简单地表示该方程组：

\[A=\begin{pmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \end {pmatrix}\]

\[B=\begin{pmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} & b_1 \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} & b_n \end {pmatrix}\]

我们将矩阵 \(A\) 称为线性方程组的系数矩阵，矩阵 \(B\) 称为线性方程组的增广矩阵。

对于最简单的情况：

\[A=\begin{pmatrix} a_{1, 1} & & & \\ & a_{2, 2} & & \\ & & \ddots & \\ & & & a_{n, n} \end{pmatrix}\]

这时只需要方程按顺序代入法求解即可。

所以要将系数矩阵转化成这样的矩阵（对角矩阵）就行了。而转化的方法就是做矩阵的初等变换（这里只需要初等行变换）：

1. 交换：\(r_i \leftrightarrow r_j\)，即交换第 \(i\) 行和第 \(j\) 行
2. 数乘：\(r_i \times k \rightarrow r_i\)，即第 \(i\) 行所有元素同时乘 \(k\)（\(k \neq 0\)）
3. 倍加：\(r_i \times k + r_j \rightarrow r_j\)，即第 \(i\) 行所有元素乘 \(k\) 加再到第 \(j\) 行

显然，对增广矩阵进行这 \(3\) 中操作并不影响原方程组的解。

（下面默认 \(r_i\) 为增广矩阵的第 \(i\) 行，对增广矩阵操作的同时对系数矩阵操作）

从左往右转化，首先找到满足 \(|a_{j,i}|\) 最大的行 \(r_j\)，交换 \(r_i\) 和 \(r_j\)。为了将系数矩阵转化成对角矩阵，第 \(i\) 列除了 \(a_{i,i}\) 全为 \(0\)，即 \(\forall j \neq i,a_{j,i}=0\)，可以用倍加操作实现，即对于所有 \(j\)，\(-r_i \times \frac {a_{i,j}} {a_{i, i}} + r_j \rightarrow r_j\)。因为在操作第 \(i\) 前已经有 \(\forall j<i, a_{i,j}=0\)，所以第 \(i\) 次操作并不会对 \(i\) 之前的列产生影响，也就是不会影响之前的操作。故 \(n\) 次操作后系数矩阵一定变成对角矩阵。

写成代码如下：

double a[MAX_N][MAX_N];					// 增广矩阵 

int main() {
	cin >> n;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n + 1; j++)
			cin >> a[i][j];						// 输入增广矩阵 
	
	for (int i = 1; i <= n; i++) {
		int maxr = i;
		for (int j = 1; j <= n; j++) 
			if (fabs(a[j][i]) > fabs(a[maxr][i]))	// 找第i列的绝对值最大的行 
				maxr = j;
		for (int j = i; j <= n + 1; j++)			// 交换操作 
			swap(a[i][j], a[maxr][j]);
		
		for (int j = 1; j <= n; j++) {
			if (j == i)
				continue; 
			double rate = a[j][i] / a[i][i]; 
			for (int k = i; k <= n + 1; k++)		// 倍加操作 
				a[j][k] -= a[i][k] * rate;
		}
	}
	
	for (int i = 1; i <= n; i++) {
		bool flag = true;
		for (int j = 1; j <= n; j++)
			flag &= fabs(a[i][j]) < eps; 
		if (flag && fabs(a[i][n + 1]) >= eps) {		// 系数全为0但和不为0，无解 
			cout << -1 << '\n';
			return 0;
		}
	}
	
	for (int i = 1; i <= n; i++)
		if (fabs(a[i][i]) < eps) {		// a[i][i]为0，x[i]绝对值最大的字数都为0，故有无穷解 
			cout << 0 << '\n';
			return 0;
		}
	
	for (int i = 1; i <= n; i++) {
		cout << fixed << setprecision(2);
		cout << a[i][n + 1] / a[i][i] << '\n';
	}
	return 0;
}

值得注意的是，有可能在转化过程中第 \(i\) 次操作时 \(\forall j>i,a_{j,i}=0\)，即根本不存在第 \(i\) 列非零的行，\(a_{i,i}\) 必然为 \(0\)，之所以取绝对值最大的行也是为了方便判断时候出现 \(a_{i,i}\) 必然为 \(0\) 的情况。但这仍有可能性。

当第 \(i\) 列全为 \(0\) 时，会有 \(maxr=i\)，正常情况下不会出问题（如果直接 continue），但在部分数据下可能影响后续操作：有可能第 \(j\) 列（\(j>i\)）中仅有 \(a_{i,j}\neq 0\)，即第 \(j\) 次操作本应该有 \(maxr=i\)，但因 \(i\) 占了本不应该占的第 \(i\) 行，导致第 \(j\) 次操作误判。

为了避免这种情况，可以直接在第 \(j\) 次操作发现存在 \(a_{i,i}=0\) 且 \(a_{i,j}>a_{maxr,j}\) 的行便直接更新 \(maxr\)，不用考虑 \(i\) 和 \(j\) 的大小关系。

改正后完整代码如下：

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-9;
const int MAX_N = 105;

int n;
double a[MAX_N][MAX_N];

int main() {
	cin >> n;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n + 1; j++)
			cin >> a[i][j];
	
	for (int i = 1; i <= n; i++) {
		int maxr = i;
		for (int j = 1; j <= n; j++)
			if (j > i || fabs(a[j][j]) < eps)
				if (fabs(a[j][i]) > fabs(a[maxr][i]))
					maxr = j;
		for (int j = i; j <= n + 1; j++)
			swap(a[i][j], a[maxr][j]);
		
		if (fabs(a[i][i]) < eps)
			continue;
		
		for (int j = 1; j <= n; j++) {
			if (j == i)
				continue; 
			double rate = a[j][i] / a[i][i];
			for (int k = i; k <= n + 1; k++)
				a[j][k] -= a[i][k] * rate;
		}
	}
	
	for (int i = 1; i <= n; i++) {
		bool flag = true;
		for (int j = 1; j <= n; j++)
			flag &= fabs(a[i][j]) < eps;
		if (flag && fabs(a[i][n + 1]) >= eps) {
			cout << -1 << '\n';
			return 0;
		}
	}
	
	for (int i = 1; i <= n; i++)
		if (fabs(a[i][i]) < eps) {
			cout << 0 << '\n';
			return 0;
		}
	
	for (int i = 1; i <= n; i++) {
		cout << fixed << setprecision(2);
		cout << a[i][n + 1] / a[i][i] << '\n';
	}
	return 0;
}

其他练习：
P3389 【模板】高斯消元法；
P5027 Barracuda；
P3164 [CQOI2014] 和谐矩阵（仅含01的异或和方程组）；
P2447 [SDOI2010] 外星千足虫（仅含01的异或和方程组）；