高斯消元——杨子曰算法

高斯消元——杨子曰算法

超链接：数学合集
（不知道这算不算数论）

高斯消元，可以干一件事情——求解n元一次方程组
黑喂狗：

我们以这样一个方程为例来讲解一下：
$$\{\begin{array}{r}2x+3y+z=14\\ 3x+y+4z=30\\ x+4y+2z=17\end{array}$$

解多元的方程其实就是不停地消元消元消元，我们有两种方法来消元：加减消元和代入消元，之后我们都会用到

我们先把上面的方程写成一个增广矩阵（←不知道是啥完全没关系，看看下面就懂了）：

$$\left[\begin{array}{cccc}2& 3& 1& 14\\ 3& 1& 4& 30\\ 1& 4& 2& 17\end{array}\right]$$

我们要注意算法的普遍适用性，毕竟我们还是要打代码的呀！

如果我们可以把这个矩阵变成这样（我们给这样的矩阵的系数部分取个名字，叫倒三角矩阵），是不是就美滋滋了：

$$\left[\begin{array}{cccc}{a}_x& a_y& a_z& s_1\\ 0& b_y& b_z& s_2\\ 0& 0& c_z& s_3\end{array}\right]\Rightarrow \{\begin{array}{r}{a}_{x}x+a_{y}y+a_{z}z=s_1\\ b_{y}y+b_{z}z=s_2\\ c_{z}z=s_3\end{array}$$
我们根据最后一个方程就可以解出最后一个z，然后把z代回中间那个方程解出y，再把y和z代入第一个方程，我们就有了x

完美解出！

根据等式的性质，我们可以对矩形进行的操作有这些：

• 用一个非零的数乘上一行
• 把其中一行的若干倍加到一行上
• 交换两行的位置

现在我们就要考虑怎样它变成倒三角矩阵捏？

我们来观察一下这个矩阵，你会发现三个等式中x只出现了一次，这就意味着我们需要选一个方程留下，把它放在第一行，并用它消掉（我们采用加减消元）其他方程的x，然后我们就可以忽视这个方程了，因为它已经被我们处理好了

如果不看这个方程的话，你会发现我们要只留下一个y，哦，那就又重复刚才的操作了，选一个方程留下y，把它放在第二行，用它消掉其他方程的y，然后再次不管这个方程

再处理z……

哦，我们的思路渐渐清晰了，现在我们要考虑一个问题，每一个未知数我们都只要留下一个方程（忽略前面我们处理好的），那么留下的是哪一个捏？理论上来说，选哪个方程留下当然不会影响结果，BUT，我们都知道电脑有一个严重的问题叫做——卡精度，为了减小精度误差，我们往往会选择当前考虑的未知数中系数绝对值最大的那一个来减小精度误差（在文章最下面我来简单讲一下为什么）

注意一点：如果某一个未知数的系数全是零，那就无解了（打代码时，我们只要看取出来的那一个绝对值最大的数是不是0即可）

我们用上面那个矩阵简单模拟一下这个过程
$$\left[\begin{array}{cccc}2& 3& 1& 14\\ 3& 1& 4& 30\\ 1& 4& 2& 17\end{array}\right]$$
我们先来解决第一个未知数x，找到系数绝对值最大的那一个，换到第一行去：
$$\left[\begin{array}{cccc}3& 1& 4& 30\\ 2& 3& 1& 14\\ 1& 4& 2& 17\end{array}\right]$$

由于我们得回要用第一行把其他行的x消去，我们就先把第一行x的系数变成1：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 2& 3& 1& 14\\ & & & \\ 1& 4& 2& 17\end{array}\right]$$

然后分别乘上其他行x的系数一一作差，把x消去：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 2-2\ast 1& 3-2\ast \frac{1}{3}& 1-2\ast \frac{4}{3}& 14-2\ast 10\\ & & & \\ 1-1\ast 1& 4-1\ast \frac{1}{3}& 2-1\ast \frac{4}{3}& 17-1\ast 10\end{array}\right]$$
我们得到了：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& \frac{7}{3}& -\frac{5}{3}& -6\\ & & & \\ 0& \frac{11}{3}& \frac{2}{3}& 7\end{array}\right]$$
我们就完美地把x消成一个了！
现在我们就可以不管第一行了
开始消y
先把绝对值最大的调上来：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& \frac{11}{3}& \frac{2}{3}& 7\\ & & & \\ 0& \frac{7}{3}& -\frac{5}{3}& -6\end{array}\right]$$
把这个式子y前的系数变成1：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& 1& \frac{2}{11}& \frac{21}{11}\\ & & & \\ 0& \frac{7}{3}& -\frac{5}{3}& -6\end{array}\right]$$
把第三个式子的y消掉：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& 1& \frac{2}{11}& \frac{21}{11}\\ & & & \\ 0& \frac{7}{3}-\frac{7}{3}\ast 1& -\frac{5}{3}-\frac{7}{3}\ast \frac{2}{11}& -6-\frac{7}{3}\ast \frac{21}{11}\end{array}\right]$$
也就是：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& 1& \frac{2}{11}& \frac{21}{11}\\ & & & \\ 0& 0& -\frac{23}{11}& \frac{345}{33}\end{array}\right]$$
然后程序还会继续执行，找到仅剩的一行把z的系数变成1：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& 1& \frac{2}{11}& \frac{21}{11}\\ & & & \\ 0& 0& 1& 5\end{array}\right]$$
于是乎我们成功地把它变成了一个到三角矩阵
我们最先得到了z=5，然后把它带回第二行，得到了：
$$\left[\begin{array}{cccc}1& \frac{1}{3}& \frac{4}{3}& 10\\ & & & \\ 0& 1& 0& 1\\ & & & \\ 0& 0& 1& 5\end{array}\right]$$
我们就有了y=1
最后我们把y和z一起带掉第一行，得到了：
$$\left[\begin{array}{cccc}1& 0& 0& 3\\ 0& 1& 0& 1\\ 0& 0& 1& 5\end{array}\right]$$
历经千辛万苦，我们得到了这个方程的解：
$$\{\begin{array}{r}x=3\\ y=1\\ z=5\end{array}$$

接下来，我来简单地解释一下为什么每次取当前未知数中系数绝对值最大的那一个式子去消其他式子中的这个未知数可以把误差降到最低

我们再假设第l个式子中未知数t的系数是$p_{lt}$

假设，我们当前在解决未知数x，它在每个式子中系数分别是$p_{1x},p_{2x},p_{3x}\cdots p_{nx}$（已经忽略考虑完的式子，下标表示第几个式子）

假设我们选择了式子i留下

我们利用上面的式子加减消元后，把其他式子下x全部消掉，第j个式子中未知数w的系数是这样算的（花点时间，结合上面的例子你一定能看懂）：
$$q_{jw}=q_{jw}-\frac{q_{jx}}{q_{ix}}\ast q_{iw}$$

于是你会发现，我们选的$q_{ix}$越大，$\frac{q_{jx}}{q_{ix}}$的期望就越小，把它看作一个单位，这个单位就越小，减完以后误差就越小！

完美！

OK，完事

c++代码（洛谷P3389）：

#include<bits/stdc++.h>
#define eps 0.0000001
using namespace std;

const int maxn=105;

int n;
double a[maxn][maxn],ans[maxn];

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]);
		}
	}
	for (int i=1;i<=n;i++){
		int r=i;
		for (int j=i+1;j<=n;j++){
			if (fabs(a[r][i])<fabs(a[j][i])) r=j;
		}			
		if (fabs(a[r][i])<eps){
			puts("No Solution");
			return 0;
		}
		swap(a[i],a[r]);
		double div=a[i][i];
		for (int j=i;j<=n+1;j++) a[i][j]/=div;
        for(int j=i+1;j<=n;j++){
            div=a[j][i];
            for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*div;
        }

	}        
	ans[n]=a[n][n+1];
	for(int i=n-1;i>=1;i--){
     	ans[i]=a[i][n+1];
        for(int j=i+1;j<=n;j++)
        ans[i]-=(a[i][j]*ans[j]);
	}
	for (int i=1;i<=n;i++){
		printf("%.2lf\n",ans[i]);
	}
	return 0;
}

参考：
https://pks-loving.blog.luogu.org/p3389-mu-ban-gao-si-xiao-yuan-fa
《算法竞赛进阶指南》 李煜东 著

于HG机房