高斯消元法详解

概念

什么是高斯消元?让我们看一看 OI-Wiki 的解释：

高斯消元法（Gauss–Jordan elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。

高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。

通俗来讲，就是解多元一次方程。

用法

高斯消元法的核心思想和它的名字一样：消元。为了方便讲解，现在让我们写出一个三元一次方程组：

$$\{\begin{array}{l}x+y+z=6\stackrel{◯}{1}\\ 5x+3y+2z=17\stackrel{◯}{2}\\ 3x+2y+z=10\stackrel{◯}{3}\end{array}$$

相信各位读者一定一眼就看得出答案，但计算机不是人，不像各位读者那么聪明（就我一个蒟蒻，呜呜呜），于是请各位读者跟随我穿越时光，回到初二的课堂。

初中老师讲过二元一次方程组的求解方法有两种：代入消元法和加减消元法。代入消元就是指用一个未知数表示出另一个未知数，然后再代入到另一个方程里，而加减消元法则是通过两个方程括倍再通过加减法消掉一个未知数。

相信各位读者已经想到高斯消元到底是什么了。好的，本期讲解完毕，我们下次再见！

开个玩笑，我们继续。

刚才我们说到二元一次方程组的求解方式有加减消元和代入消元，那我们现在把它们扩展到 $n$ 元一次方程组，结果会怎么样呢？

答案：代入消元太复杂，加减消元得高消。

这里的高消就是指高斯消元。

很明显，代入消元对计算机和程序员来讲难度太大了，因此，一位伟大的数学家——高斯，便发明了高斯消元（不知他的初衷是啥）！

高斯消元的做法是这样的：首先把每一项的系数和右边的常数项提取出来构成矩阵，比如上面的那个方程就可以被转化成这样：

$$\{\begin{array}{l}x+y+z=6\\ 5x+3y+2z=17\\ 3x+2y+z=10\end{array}\to \left[\begin{array}{cccc}1& 1& 1& 6\\ 5& 3& 2& 17\\ 3& 2& 1& 10\end{array}\right]$$

接着让我们把第 $i$ 行第 $i$ 列的数变成 $1$，也就是给第 $i$ 行的每个数除以 $a_{i,i}$，在上面的矩阵中我们先处理第一行：

$$\left[\begin{array}{cccc}1& 1& 1& 6\\ 5& 3& 2& 17\\ 3& 2& 1& 10\end{array}\right]\to \left[\begin{array}{cccc}1& 1& 1& 6\\ 5& 3& 2& 17\\ 3& 2& 1& 10\end{array}\right]$$

然后我们以第 $i$ 行第 $i$ 列的数为基准，把第 $i$ 列第 $i$ 行下面的所有数变成 $0$，这一步解释起来有点复杂，大家自己看图：

$$\begin{array}{rl}& \left[\begin{array}{cccc}1& 1& 1& 6\\ 5& 3& 2& 17\\ 3& 2& 1& 10\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 5-5\times 1& 3-5\times 1& 2-5\times 1& 17-5\times 6\\ 3-3\times 1& 2-3\times 1& 1-3\times 1& 10-3\times 6\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& -2& -3& -13\\ 0& -1& -2& -8\end{array}\right]\end{array}$$

最后重复上述操作，直到形成一个上三角结构，为了方便读者理解，这里把全过程展示了出来：

$$\begin{array}{rl}& \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& -2& -3& -13\\ 0& -1& -2& -8\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& -2÷(-2)& -3÷(-2)& -13÷(-2)\\ 0& -1& -2& -8\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& -1& -2& -8\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& -1-1\times (-1)& -2-\frac{3}{2}\times (-1)& -8-\frac{13}{2}\times (-1)\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& 0& -\frac{1}{2}& -\frac{3}{2}\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& 0& -\frac{1}{2}÷(-\frac{1}{2})& -\frac{3}{2}÷(-\frac{1}{2})\end{array}\right]\\ & \Rightarrow \left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& 0& 1& 3\end{array}\right]\end{array}$$

于是我们就得到了一个上三角结构的矩阵（即下半部分的三角形全是 $0$），而它的等效方程就是这个：

$$\left[\begin{array}{cccc}1& 1& 1& 6\\ 0& 1& \frac{3}{2}& \frac{13}{2}\\ 0& 0& 1& 3\end{array}\right]\iff \{\begin{array}{l}x+y+z=6\\ y+\frac{3}{2}z=\frac{13}{2}\\ z=3\end{array}$$

接下来就是最重要的一步了：代入消元。有人会问：坐着你前面不是说代入消元很麻烦吗，怎么现在又要代入消元了？是不是玩不起啊？

其实不是的，因为通过观察左边的那个方程，我们会发现 $z$ 的值实际上已经解开了，对应到矩阵上就是最后一行的倒数第二个是一个 $1$，而前面的数全是 $0$，这时最后一行的最后一个数就是最后的那个未知数的解，记录下来。

然后回到倒数第二行，我们不难发现：除了倒数第二个数（注意：我这里没写第二个，而写的是倒数第二个）是 $1$ 之外，其他的不是 $0$ 就是其他乱七八糟的数，在右边的方程中，我们在把 $z$ 解开后肯定是带入到倒数第二个方程中，解开 $y$，而对应在左边的矩阵里就是用最后一个数减去最后那个数乘以下面那个未知数的值，于是得到了 $y=2$ 这个解，记录下来。

然后就不需要我讲了吧。（不会有人不会举一反三吧。）

特殊情况

上面的操作中，由于我们过分理想化的假设，导致我们没遇到一种特殊情况：系数为 $0$！我们举个例子看看会发生什么：

$$\{\begin{array}{l}y+z=3\\ x+z=4\\ x+y=5\end{array}$$

画成矩阵：

$$\left[\begin{array}{cccc}0& 1& 1& 3\\ 1& 0& 1& 4\\ 1& 1& 0& 5\end{array}\right]$$

然后按照上面的步骤：首先我们需要将 $a_{1,1}$ 变成 $1$。但我们发现个事：$a_{1,1}=0$！按照正常的做法只需要除以一个 $a_{1,1}$ 就行了，但是没有除以 $0$ 这个东西啊喂！于是你的程序会直接 $RE$ 掉。

但是这个方程没有解吗？并不是。上过小学奥数或者初中的人都知道这个方程有解 $\{\begin{array}{l}x=3\\ y=2\\ z=1\end{array}$（呀，一不小心说出答案了），那这是怎么一回事呢？

原因很简单：只是因为这个方程刚好系数为 $0$。你看下一个方程系数不就又是 $1$ 了吗？因此我们可以适当把第二行和第一行的数换个位置，于是这个矩阵变成了这样：

$$\left[\begin{array}{cccc}1& 0& 1& 4\\ 0& 1& 1& 3\\ 1& 1& 0& 5\end{array}\right]$$

然后就可以快乐的解方程了。

总结一下： 规则就两条：

1. 能直接解的直接解。
2. 不能直接解的换两行再解。

最后说明一下判无解和无数解的办法：无解就是存在某一行除了最后一列其他全是 $0$，无数解就是存在一行全是 $0$。

浅浅算一下时间复杂度：$O(n^3)$。

代码太简单，就不展示了。（我绝对不会告诉你是因为我不会写。）

但是！

---

我的奇怪方法

由于本蒟蒻太菜，实在不会写高消，于是我自己再高消的基础上改了一下。

首先，我不需要换行，对于每一行我都遍历一遍找一个主元，而不是一定以第 $i$ 行的第 $i$ 个为主元。

其次，我每次消是把这一列全消成 $0$（除这个主元之外），这样最后就是一个 $01$ 矩阵。

最后，我不需要还原，因为左边是个 $01$ 矩阵，所以系数只有可能是 $0$ 或 $1$，所以直接输出最后一列就行。

当然，我的解是随机分布的，因此需要拿个标记数组标记每个解再第几行再输出。

时间复杂度：$O(n^3)$。

代码：

//洛谷P3389
#include<bits/stdc++.h>
#define int long long
#define code using
#define by namespace
#define plh std
code by plh;
int n,vis[106],h[106];
const double eps=1e-7;
double a[106][106],f[106];
void Gauss()
{
	for(int i=1;i<=n;i++)
	{
		int id=0;
		for(int j=1;j<=n;j++)
		{
			if(fabs(a[i][j])>eps&&!vis[j])
			{
				id=j;
				break;
			}
		}
		if(id)
		{
			vis[id]=1;
			h[id]=i;
			for(int j=1;j<=n;j++)
			{
				if(j!=id)
				{
					a[i][j]/=a[i][id];
				}
			}
			f[i]/=a[i][id];
			a[i][id]=1;
			for(int j=1;j<=n;j++)
			{
				if(j!=i)
				{
					for(int k=1;k<=n;k++)
					{
						if(k!=id)
						{
							a[j][k]-=a[i][k]*a[j][id];
						}
					}
					f[j]-=f[i]*a[j][id];
					a[j][id]=0;
				}
			}
		}
		else
		{
			if(fabs(f[i])>eps)
			{
				cout<<"No Solution";
				n=0;
				return;
			}
		}
	}
}
signed main()
{
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			cin>>a[i][j];
		}
		cin>>f[i];
	}
	Gauss();
	if(!n)
	{
		return 0;
	}
	for(int i=1;i<=n;i++)
	{
		if(!vis[i])
		{
			cout<<"No Solution";
			return 0;
		}
	}
	for(int i=1;i<=n;i++)
	{
		printf("%0.2lf\n",f[h[i]]);
	}
	return 0;
}

说句实话，这篇文章制作非常不易（特别是 LateX 那一部分，真的是一个字一个字打出来的），希望读者能点个赞，非常感谢！