p2455 线性方程组（高斯消元）

p2455才是高斯消元的真模板题好吧

高斯消元

题面：传送门

对不起我真的是太懒了

先来明确高斯消元是干什么的：

高斯消元是一种求解线性方程组的方法。

所谓线性方程组，就是由\(M\)个\(N\)元一次方程组共同构成的。

线性方程组的所有系数可以写成一个\(M\)行\(N\)列的系数矩阵，再加上每个方程等号右侧的常数，可以写成一个\(M\)行\(N+1\)列的“增广矩阵”。

\(\begin{cases}x_1+2x_2-x_3=-6\\2x_1+x_2-3x_3=-9\\-x_1-x_2+2x_3=7\end{cases}\Rightarrow\left[\begin{array}{ccc|c}1 &2&-1&-6\\2&1&-3&-9\\-1&-1&2&7 \end{array}\right] \)

相信各位小学或者初中一定学过二元一次方程组吧

如果没学过那我觉得宁可以退役了

有个消元的方法相信你记忆犹新：加减消元法；

其实高斯消元法干的事情和加减消元法是差不多的；

步骤大概就是：

我们将一行加上另一行的若干倍，这样就可以消去其中一行的某一个系数；

我们可以来看一下下面的过程：

\(\begin{bmatrix} 1&2&-1&-6\\2&1&-3&-9\\-1&-1&2&7\end{bmatrix}\)

我们将第二行每一项分别加上第一行的\(-2\)倍，于是我们发现，第二行的第一项系数被消成\(0\)了。如下

\(\begin{bmatrix} 1&2&-1&-6\\0&-3&-1&3\\-1&-1&2&7\end{bmatrix}\)

接下来我们都这么操作。

对于第\(n\)行，我们让该行的第\(n\)个元素作为“主元”，并且将其他行的同列的系数都消掉，最后我们可以得到下面的形式：

\(\left[ \begin{array}{ccc|c} 1&0&0&1\\0&1&0&-2\\0&0&1&3 \end{array} \right]\)

是的，成了这样以后，发现主对角线上的数和常数一般都不是\(0\)，其他位置都是\(0\)了（后面还有一些特殊情况与判定技巧），我们就可以愉快的写出答案：

\(x_1=1\)，\(x_2=-2\)，\(x_3=3\)；

是的，这就是高斯消元算法。

但是在实际应用中，我们也许会碰到很多奇奇怪怪的情况。

就例如你只掌握上面说到的东西是远远不够通过这道题的。

题目中还有一些更高端的东西：什么时候有无穷多的实数解？什么时候无解？这都是我们需要考虑的；

首先，在高斯消元的过程中，可能出现\(0=d\)这样的方程，其中\(d\)是一个非零常数，这表明这些方程中出现了矛盾，方程组无解。

其次，我们也可能出现这样的情况：

我们高斯消元结束以后，发现某一行系数全都是\(0\)，并且常数也是\(0\).

啊就是出现了\(0==0\)的情况？？

是的，此时方程组有无数解，如下图：

\(\left[ \begin{array}{ccc|c} 1&2&0&4\\0&0&1&1\\0&0&0&0 \end{array} \right]\)

那么我们总结一下：

1.当高斯消元结束以后，如果存在系数全是\(0\)，但是常数不是\(0\)的行，那么方程组无解

2.在高斯消元结束以后，如果存在系数全是\(0\)，并且常数也是\(0\)的行，那么方程组有无数解

学会了这些就能通过这些题了吗？

当然不是！

还有一个东西你没有注意到：精度损失

精度损失是怎么来的?如何降低精度损失呢?

我们在加减消元的时候，有的时候需要加的也许不是一个整数，可能是某一行的几分之几倍。

于是，考虑到这，我们发现精度损失的原因了。

所以，我们就要让产生几分之几倍（确切的说：当式子一的系数是\(a_1\)，式子二的系数是\(a_2\)时，我们加的这个数就是\(a_2/a_1\)倍（注意我们主元系数是在分母上））的这个过程有较小的精度损失

我们可以思考一下，我们在消元的时候一定要刻意按照题目中给的顺序一行一行来嘛？

也许可以不这样。

在什么情况下，我们能产生较小的精度损失呢？

我们举个例子来看：

考虑这样的一个方程组：

\(1000x_1+ax_2+bx_3=c\)）

\(0.5x_1+dx_2+ex_3=f\)

我们设精度误差\(eps=1e-8\)；

如果我们选择\(1000\)作为主元，

那么理论值\(ima_1=0.5/1000\)

实际值\(rea_1=0.5+eps/1000+eps\)

误差值Δ在\(rea_1-ima_1=\)一个超级小的数

---

如果我们选择\(0.5\)作为主元，那么理论值\(ima_2=1000/0.5\)

实际值\(rea_2=1000+eps/0.5+eps\)；

误差值Δ在\(rea_2-ima_2=40\)左右

发现差距了吗？？

知道该怎么做了吗？？

是的，我们每次选择较大的作为主元就好了！

代码实现嘛，长这样的;

	double maxx;//最大值
	int maxline;//最大值所在行
	for(re int i=1;i<=n;++i){//枚举每一列
		maxx=0,maxline=i;
		for(re int j=1;j<=n;++j)//枚举每一行的同一列位置
			if(!vis[j]&&abs(a[j][i])>maxx){//如果没来过并且很大（因为我们加的数正负都可以，所以我们取绝对值
				maxx=abs(a[j][i]),maxline=j;//记录最大值以及所在的行数
			}
            
		//........
        
		if(maxline!=i)swap(a[i],a[maxline]);//如果发现我们在消元过程中，这个位置不是需要的最大的，我们就把这一行和最大的那一行换一下位置
		vis[i]=1;//来过了
	}	

掌握了这些，这个题就可以\(A\)了

CODE：

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

#define re register

const int N=105;
const double eps=1e-8;//精度误差

double a[N][N];//矩阵
bool vis[N];
int n;

inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}
	return x*f;
}

int main(){
	n=read();
	
	for(re int i=1;i<=n;++i)
		for(re int j=1;j<=n+1;++j)a[i][j]=read();//读入矩阵，把常数也读进去所以是n+1
		
	double maxx;
	int maxline;
	for(re int i=1;i<=n;++i){
		maxx=0,maxline=i;
		for(re int j=1;j<=n;++j)
			if(!vis[j]&&abs(a[j][i])>maxx){//找最大的
				maxx=abs(a[j][i]),maxline=j;
			}
		
		if(maxx<=eps)continue;/如果太小了，我们可以认为是0，并且把他跳过去，因为已经是0了不需要再操作
		if(maxline!=i)swap(a[i],a[maxline]);//将系数最大的换到i行
		vis[i]=1;
		
		for(re int j=1;j<=n;++j){//对其他行关于i列进行消元
			if(j!=i){
				double t=a[j][i]/a[i][i];
				for(re int k=i;k<=n+1;++k){
					a[j][k]-=t*a[i][k];
				}
			}
		}
	} 
	
	bool flag=0;//因为在消元完成后，只有主对角线上系数不是零，所以我们只需要检查一下，如果主对角线上出现了0的那一行的常数是不是0就可以了
	for(re int i=1;i<=n;++i){
		if(abs(a[i][i])<=eps){
			flag=1;//如果发现主对角线上有0存在，那么一定出了特殊情况，要么是无解，要么是无穷解
			if(abs(a[i][n+1])>eps){//如果常数不是零，那么可以认定就是无解
				puts("-1");
				return 0;
			}
		}
	}
		
	if(flag){
		puts("0");//如果主对角线上有0，并且刚才检测过，常数肯定是0，那么就一定有无穷多解
		return 0;
	}
	
	for(re int i=1;i<=n;i++){//不然我们就要输出每一项的解是多少
		printf("x%d=%.2lf\n",i,a[i][n+1]/a[i][i]);//这里解释一下，为什么要用常数除以系数：因为我们在消元以后，系数并不一定变成了1，也可能是其他的常数，那么就成了kx=b的形式，所以我们拿b/k就是x了
	}
	
	return 0;
}