【数学】高斯-约旦消元法

给定 $n$ 元一次方程组
$$\{\begin{array}{l}{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{array}$$
请求出方程组的解的情况：

• 无解；

• 无穷多解；

• 唯一解。

对于这样的问题，我们可以使用 高斯消元法 进行求解，当然高斯消元法有一个回代的过程，代码略长，而且精度较低。

所以我们隆重推出 高斯-约旦消元法 ！！！

回顾一下我们是怎么手算的，一般用的都是 加减消元法，普通高斯和高斯-约旦用的都是加减消元。

在此之前，我们需要了解一下矩阵初等变换。

在线性代数中，矩阵初等行变换 是指以下三种变换类型 ：

1. 交换矩阵的两行；

2. 用一个非零数 $k$ 乘矩阵的某一行所有元素；

3. 把矩阵的某一行所有元素乘以一个数 $k$ 后加到另一行对应的同一列的元素上；

类似地，把以上的 行 改为 列 便得到 矩阵初等列变换 的定义。

矩阵初等行变换与初等列变换合称为 矩阵初等变换。

若矩阵 $A$ 经过有限次的初等行变换变为矩阵 $B$，则矩阵 $A$ 与矩阵 $B$ 行等价；若矩阵 $A$ 经过有限次的初等列变换变为矩阵 $B$，则矩阵 $A$ 与矩阵 $B$ 列等价；若矩阵 $A$ 经过有限次的初等变换变为矩阵 $B$，则矩阵 $A$ 与矩阵 $B$ 等价。

当然列的用不着

首先有一个由系数构成的 $n\times n$ 的矩阵
$$\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}\end{array}\right]$$

然后是一个由常数构成的 $n\times 1$ 的列向量
$$\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

把它们放在一起构成一个 $n\times (n+1)$ 的增广矩阵：
$$\left[\begin{array}{cccccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}& \mid & b_1\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}& \mid & b_2\\ ⋮& ⋮& \ddots & ⋮& \mid & ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}& \mid & b_n\end{array}\right]$$
我们遍历每一列，对于第 $i$ 列选出最大的 未处理过的行 上的数作为主元，意味着加减消元后除了主元这一行外其他行的第 $i$ 列都为 $0$。

选最大的作为主元的原因是避免精度误差。

然后就是加减消元了。

举个例子
$$\{\begin{array}{l}2x-y+z=1\\ 4x+y-z=5\\ x+y+z=0\end{array}$$
增广矩阵就是
$$\left[\begin{array}{ccccc}2& -1& 1& \mid & 1\\ 4& 1& -1& \mid & 5\\ 1& 1& 1& \mid & 0\end{array}\right]$$
先是第 $1$ 列，选 $4$ 为主元。

$4$ 在第 $2$ 行，将第 $2$ 行与第 $1$ 行交换。
$$\left[\begin{array}{ccccc}4& 1& -1& \mid & 5\\ 2& -1& 1& \mid & 1\\ 1& 1& 1& \mid & 0\end{array}\right]$$

对于第 $2$ 行，第一列上的 $2$ 是 $4$ 的 $\frac{1}{2}$。
$$\left[\begin{array}{ccccc}4& 1& -1& \mid & 5\\ 2-4\times \frac{1}{2}& -1-1\times \frac{1}{2}& 1-(-1)\times \frac{1}{2}& \mid & 1-5\times \frac{1}{2}\\ 1& 1& 1& \mid & 0\end{array}\right]$$
化简得
$$\left[\begin{array}{ccccc}4& 1& -1& \mid & 5\\ 0& -\frac{3}{2}& \frac{3}{2}& \mid & -\frac{3}{2}\\ 1& 1& 1& \mid & 0\end{array}\right]$$
对第 $3$ 行的处理同理
$$\left[\begin{array}{ccccc}4& 1& -1& \mid & 5\\ 0& -\frac{3}{2}& \frac{3}{2}& \mid & -\frac{3}{2}\\ 0& \frac{3}{4}& \frac{5}{4}& \mid & -\frac{5}{4}\end{array}\right]$$

现在到了第 $2$ 列，未处理过的是 $2,3$ 行，选最大的 $\frac{3}{4}$ 为主元。

将第 $3$ 行与第 $2$ 行交换
$$\left[\begin{array}{ccccc}4& 1& -1& \mid & 5\\ 0& \frac{3}{4}& \frac{5}{4}& \mid & -\frac{5}{4}\\ 0& -\frac{3}{2}& \frac{3}{2}& \mid & -\frac{3}{2}\end{array}\right]$$

消元得
$$\left[\begin{array}{ccccc}4& 0& -\frac{8}{3}& \mid & \frac{20}{3}\\ 0& \frac{3}{4}& \frac{5}{4}& \mid & -\frac{5}{4}\\ 0& 0& 4& \mid & -4\end{array}\right]$$
第 $3$ 列，未处理过的是第 $3$ 行，选 $4$ 作主元。
$$\left[\begin{array}{ccccc}4& 0& 0& \mid & 4\\ 0& \frac{3}{4}& 0& \mid & 0\\ 0& 0& 4& \mid & -4\end{array}\right]$$

这样就把原来的矩阵通过初等变换，使得系数矩阵只有主对角线位置的元素非零，即一个对角矩阵。

上面那个矩阵的意思是
$$\{\begin{array}{l}4x=4\\ \frac{3}{4}y=0\\ 4z=-4\end{array}$$
所以再把系数除过去就得到
$$\{\begin{array}{l}x=1\\ y=0\\ z=-1\end{array}$$
请自行检验。

时间复杂度为 $\mathcal{O}(n^3)$。

当然方程还可能是无解或无穷多解。

考虑一个一元一次方程 $ax=b$ 的解的情况：

1. 无解：$a=0,b\ne 0$；
2. 无穷多解：$a=b=0$；
3. 唯一解：$a\ne 0$。
4. 所以当发现某一列的主元为 $0$ 时，因为主元是最大的，所以意味着这一列全都是 $0$，那么要么无解，要么有无穷多解。

P3389 【模板】高斯消元法

$\text{Code}$

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 105;

int n;
db a[MAXN][MAXN];

bool Gauss_Jordan()
{
	for (int i = 1; i <= n; i++)
	{
		int mx = i;
		for (int j = i + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) < fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != i)
		{
			swap(a[i], a[mx]);
		}
		if (!a[i][i])
		{
			return false;
		}
		for (int j = 1; j <= n; j++)
		{
			if (i != j)
			{
				db val = a[j][i] / a[i][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[i][k] * val;
				}
			}
		}
	}
	return true;
}

int main()
{
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n + 1; j++)
		{
			scanf("%lf", a[i] + j);
		}
	}
	if (Gauss_Jordan())
	{
		for (int i = 1; i <= n; i++)
		{
			printf("%.2lf\n", a[i][n + 1] / a[i][i]);
		}
	}
	else
	{
		puts("No Solution");
	}
	return 0;
}

关于判断无解和无穷多解

P2455 [SDOI2006]线性方程组

这下要具体到到底是无解还是无穷多解了。

这就是高斯-约旦的一个缺点：相比于普通高斯，它更难判断无解和无穷多解。

但是也是可以处理的。

具体地，我们用 $r$ 来记录当前行，如果主元为 $0$，那么 $\mathrm{continue}$ 到下一列，但 $r$ 不变；否则消元后令 $r\leftarrow r+1$。

遍历完所有列后：

• 若 $r=n$，说明有唯一解；
• 若 $r<n$，说明第 $r+1\sim n$ 行系数矩阵全都是 $0$，若列向量全是 $0$，说明有无穷多解，否则无解。

$\text{Code}$

#include <iostream>
#include <cstdio>
#include <cmath>
typedef double db;
using namespace std;

const int MAXN = 55;

int n;
db a[MAXN][MAXN], ans[MAXN];

int Gauss_Jordan()
{
	int r = 1;
	for (int i = 1; i <= n; i++)
	{
		int mx = r;
		for (int j = r + 1; j <= n; j++)
		{
			if (fabs(a[j][i]) > fabs(a[mx][i]))
			{
				mx = j;
			}
		}
		if (mx != r)
		{
			swap(a[r], a[mx]);
		}
		if (!a[r][i])
		{
			continue;
		}
		for (int j = 1; j <= n; j++)
		{
			if (j != r)
			{
				db val = a[j][i] / a[r][i];
				for (int k = i + 1; k <= n + 1; k++)
				{
					a[j][k] -= a[r][k] * val;
				}
			}
		}
		r++;
	}
	if (--r < n)
	{
		for (int i = r + 1; i <= n; i++)
		{
			if (a[i][n + 1])
			{
				return -1;
			}
		}
		return 0;
	}
	for (int i = 1; i <= n; i++)
	{
		ans[i] = a[i][n + 1] / a[i][i];
		if (!ans[i])
		{
			ans[i] = 0;
		}
	}
	return 1;
}

int main()
{
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n + 1; j++)
		{
			scanf("%lf", a[i] + j);
		}
	}
	int res = Gauss_Jordan();
	if (res != 1)
	{
		printf("%d\n", res);
	}
	else
	{
		for (int i = 1; i <= n; i++)
		{
			printf("x%d=%.2lf\n", i, ans[i]);
		}
	}
	return 0;
}