[算法] 高斯消元详解

原文链接（强烈建议安利）

0.前置知识

• 知道如何解三元一次方程组
• 有手，有脑子

1.答案的表示与存储

先解一个方程组：

2x+3y+5z=31
x-4y -z=-6
4x+2y-5z=9

---

我们把这个方程组写成 机器能读懂 的 表格形式 ：

2 3 5 31
1 -4 -1 -6
4 2 -5 9

第一列代表 $x$ 的系数，第二列代表 $y$ 的系数……
注意多出来的第四列是答案的具体数值。
第一行代表式子 1，第二列代表式子 2……

---

因为老是使用 $x,y,z$ 表示 $n$ 个数不方便，所以我们使用 $a_1,a_2...a_n$ 来表示。

即，上面的表格的意义就是：

第 $i(1\le i\le n)$ 列代表每个式子中 $a_i$ 的系数，第 $n+1$ 列表示每个式子的具体数值。（因为存储的时候这个都是普通的数值，所以读者务必辨清）

第 $i(1\le i\le n)$ 行代表第 $i$ 个式子。

接着全文中，使用 $a[i][j]$ 表示第 $i$ 行第 $j$ 列的数值/系数。

---

2.高斯消元思想

名字这么高级，但是方法很土。具体操作是：加减消元法+代入消元法。

首先，接着上面的系数矩阵：

2 3 5 31
1 -4 -1 -6
4 2 -5 9

然后开始高斯消元。

枚举 $i$ ，顺序消除未知数 $a_i$ 。

然后枚举 $j$ ，找到一个绝对值最小不为 $0$ 的 $a[j][i]$ 。

找到那个构成最小 $a[j][i]$ 的 $j$ 后，将第 $j$ 个式子与第 $i$ 个式子交换

枚举 $j$ ，设 $mul=a[i][i]/a[j][i]$ ，
将 $j$ 式子乘上 $mul$ 并且减去 $i$ 式子

我们一层层剖析为什么这样做。
首先，枚举 $i$ ，需要消掉第 $i$ 个未知数。
然后，选择 $j$ 使得 $a[j][i]$ 最小，这个其实不是必须的（后面讲为什么）。交换 $j$ 式子和 $i$ 式子。
再后来，枚举每一个 $j$ ，将 $j$ 式子乘上 $a[i][i]/a[j][i]$ ，是因为需要将 $j$ 式子第 $i$ 项的系数 调整与 $i$ 式子第 $i$ 项的系数一样；
将 $j$ 式子减去 $i$ 式子，目的为了消去第 $i$ 项的系数 。

为什么要选择 绝对值最小而不为 $0$ 的 $a[j][i]$ 的 $j$ 呢？因为绝对值最小的 $a[j][i]$ 大致可以构成更多整除的情况，而这样就可以避免掉精度。
不能取 $0$ 是当然的。

如果还没有理解？那么一层层模拟。
还是以这个方程组为例子：

2 3 5 31
1 -4 -1 -6
4 2 -5 9

首先，消去第一个系数。

1 -4 -1 -6
2 3 5 31
4 2 -5 9

找到最小 $a[1]$ 的式子，调到第一个；

1 -4 -1 -6
1 3/2 5/2 31/2
4 2 -5 9

第二个柿子 $\times 0.5$ 。

1 -4 -1 -6
0 11/2 7/2 43/2
4 2 -5 9

第二个式子减去第一个。

---

以此类推，最后这个答案的系数矩阵是：（没有对齐，有点丑）

4 0 0 20
0 -4.5 0 -9
0 0 7.611 22.833

然后最后的 $n+1$ 系数除以 第 $i$ 行第 $i$ 列的数值即可。

3.你们想看的代码。

P3389 【模板】高斯消元法

注释有点少，能看懂就行。

#include<bits/stdc++.h>
#define RI register int
using namespace std;typedef long long LL;const int inf=1073741823;int FL,CH;template<typename T>void in(T&a){for(a=0,FL=1,CH=getchar();!isdigit(CH);CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH);CH=getchar())a=a*10+CH-'0';a*=FL;}template<typename T,typename ...Args>void in(T&a,Args&...args){in(a);in(args...);}
const int maxn=110;const double eps=1e-6;
double a[maxn][maxn];
int n;
int main()
{
	cin>>n;
	for(RI i=1;i<=n;++i){
		for(RI j=1;j<=n+1;++j)scanf("%lf",&a[i][j]);}
	//输入系数矩阵
	for(RI rp,i=1;i<=n;++i)
	{
		rp=i;
		for(RI j=i+1;j<=n;++j)
		if(fabs(a[j][i])>fabs(a[rp][i]))rp=j;
		if(fabs(a[rp][i])<=eps)
		return printf("No Solution"),0;
		//寻找绝对值最小的 j
		//这个特殊一点，要判无解
		if(i!=rp)swap(a[rp],a[i]);
		//交换 j 和 i
		double div=a[i][i];
		for(RI j=1;j<=n;++j)
		if(j!=i){
			div=a[j][i]/a[i][i];
			//计算 mul 值
			for(RI k=1;k<=n+1;++k)
				a[j][k]-=a[i][k]*div;
			//乘上，再减去
			//高斯消元算法的精髓
		}
	}
	for(RI i=1;i<=n;++i)
		printf("%.2lf\n",a[i][n+1]/a[i][i]);//回代，并输出答案
	return 0;
}