2021-07-19

高斯消元

*原题：题目背景
输入一个包含 n 个方程 n 个未知数的线性方程组。

方程组中的系数为实数。

求解这个方程组。

下图为一个包含 m 个方程 n 个未知数的线性方程组示例：

输入格式
第一行包含整数 n。

接下来 n 行，每行包含 n+1 个实数，表示一个方程的 n 个系数以及等号右侧的常数。

输出格式
如果给定线性方程组存在唯一解，则输出共 n 行，其中第 i 行输出第 i 个未知数的解，结果保留两位小数。

如果给定线性方程组存在无数解，则输出 Infinite group solutions。

如果给定线性方程组无解，则输出 No solution。

数据范围
1≤n≤100,
所有输入系数以及常数均保留两位小数，绝对值均不超过 100。

输入样例：
3
1.00 2.00 -1.00 -6.00
2.00 1.00 -3.00 -9.00
-1.00 -1.00 2.00 7.00
输出样例：
1.00
-2.00
3.00

解题思路：利用高斯消元法，将系数矩阵转换为阶梯矩阵，时间复杂度为O(n³)。阶梯矩阵就是每行的第一个非0元素为1，或者该行全部为0（系数矩阵该行全部为0）。为方便计算，需要将答案作为单独一列加入到系数矩阵

转换为阶梯行列式的步骤
1-> 为方便计算，找到某一列绝对值最大的元素的那一个行，将该行作为操作行，将其置顶，如果操作行数为i，那么需要将该行第一个非0元素转换为1，并将该元素下面的所有元素清0，通过方程式变形和两两之间的运算来修改未知数的系数（详情可参考线性代数矩阵变为阶梯矩阵）

以下为操作代码

	int c,r;  //表示行和列 
	//化解行列式，将行列式化为最简阶梯行列式 
	for(c=0,r=0;c<n;c++)   //遍历每一列，将该列绝对值最大的数所在行置顶，方便计算
	{
		int t=r;  //记录最大绝对值元素所在行 
		for(int i=r;i<n;i++)  //每次操作还未确立的行,找到此列绝对值最大的元素所在行，将其记录 
		    if(fabs(g[i][c])>fabs(g[t][c]))  t=i;  //寻找最大的数值是因为可以避免系数变得太大，精度较高.
		    
		if(fabs(g[t][c])<eps) continue; //绝对值最大等于(接近)0，那么所有该列数都是0,无需进一步操作
		for(int i=c;i<n+1;i++)   //否则从第c列开始交换行，从第c列开始是操作完成的数据就无需再操作了
		swap(g[t][i],g[r][i]);   //r为最顶行，t为最大绝对值所在行，按列交换
		for(int i=n;i>=c;i--)   //将最大绝对值元素所在行转换为最简行列式(阶梯行列式)，第一个非0元素，变为1，其余元素除以相同的值 
		g[r][i]/=g[r][c]; //到过来算是因为没有存储使其变为1的约数，顺序算会导致本身变为1，则后面的元素都会除以1 
		for(int i=r+1;i<n;i++)  //在把该列该行下面的所有元素消为0
		    if(fabs(g[r][c])>eps)  //该元素不为0，则操作，否则无需操作
			for(int j=n;j>=c;j--)  //方程两边进行相同的操作 
			g[i][j]-=g[r][j]*g[i][c];  //每次操作是减掉此行第一个非0元素与上一行第一个非0元素的积，这样就能保住此行第一个非0
			                           //元素变为0 
	        r++;   //操作完一轮之后，r++，原来的最上面一行就操作完成了 
	}

当我们推出 r<n时，矩阵则有无解和无穷多解
因为r<n代表有一行系数全为0，那么等式左边为0，如果等式右边为0，那么不管未知数怎么取，等式都成立，即有无穷多解。若右边不等于0，则等式矛盾，不管未知数怎么取值，等式都不成立，即无解

	//所有行和列转换完成后，判断是否有解
	if(r<n)  //说明剩下的方程个数是小于n的，解不为一，判断是无解还是无穷多解
    {
    	for(int i=r;i<n;i++)  //r及一下的面的行，由于已经转换为阶梯行列式，所以所有系数都为0,即等式左边全为0，若等式右边
		{                      //不为0，那么左边不等于右边，则无解，否则左边等于右边，有无穷多组解
		    if(fabs(g[i][n])>eps) 
		    return 0;  //无解 
		    else 
			return 2;	//无穷多解 
		}
	        
	}

如果r=n,代表n个未知数有n个一次方程，那么有唯一解
关于如何求解每个未知数，我们只需消掉后面的未知数即可，如果是第1行，则消掉第二个到最后一个元素，第二行则消掉第3个到最后一个元素，依次类推，且从倒数第二行向上推，最后一个无需推理，在先前计算阶梯矩阵时，答案已经被计算出来。
每一行为了消去应该消去的未知数，建立循环，每次消去下一行和下一个元素

	else   //r==n  有唯一解，从下往上推出每个未知数
	{
		for(int i=n-1;i>=0;i--)  //从下往上遍历每一行,当前这一行减去下一行与当前这一行的相应系数，便于消除一个未知数 
		   for(int j=i+1;j<n+1;j++)
		   g[i][n]-=g[j][n]*g[i][j];  //等式左右两边进行相同操作
		                              //如两个方程 x1+3*x2=d1  x2=d2,那么 将1式减去2式*3则可以得出x1答案，因为是阶梯行列式
									  //所以需要依次向上推，最后一行无需推，最后一行由于先前计算阶梯矩阵已经计算出结果了 
		return 1;  //唯一解 
    }    
	      

每次为了消去相应元素，如第一行消元，第一次循环应该是消去第2个元素，即行数为i，减j=i+1行第j=i+1个元素，直到循环到j=n，而i行的答案矩阵的元素减去下一行的答案矩阵的元素和当前行第j个元素的成绩（如x1+3x2=d1,x2=d2,那么x1=d2-d13）（自行根据代码推导更容易理解）

完整代码

#include <iostream>
#include <cmath>
#include <algorithm> 
using namespace std;
const int N=200;
double g[N][N]; //系数矩阵，从00开始，g[0][0]表示第一行x1的系数 
int n;
const double eps = 1e-6;  //误差不超过eps，当a<eps时，可以将a认为为0,不直接写等于0是防止计算机在运算时精度丢失 
int gauss()  //返回0表示无解，返回1表示唯一解，返回2表示无穷多解 
{
	int c,r;  //表示行和列 
	//化解行列式，将行列式化为最简阶梯行列式 
	for(c=0,r=0;c<n;c++)   //遍历每一列，将该列绝对值最大的数所在行置顶，方便计算
	{
		int t=r;  //记录最大绝对值元素所在行 
		for(int i=r;i<n;i++)  //每次操作还未确立的行,找到此列绝对值最大的元素所在行，将其记录 
		    if(fabs(g[i][c])>fabs(g[t][c]))  t=i;  //寻找最大的数值是因为可以避免系数变得太大，精度较高.
		    
		if(fabs(g[t][c])<eps) continue; //绝对值最大等于(接近)0，那么所有该列数都是0,无需进一步操作
		for(int i=c;i<n+1;i++)   //否则从第c列开始交换行，从第c列开始是操作完成的数据就无需再操作了
		swap(g[t][i],g[r][i]);   //r为最顶行，t为最大绝对值所在行，按列交换
		for(int i=n;i>=c;i--)   //将最大绝对值元素所在行转换为最简行列式(阶梯行列式)，第一个非0元素，变为1，其余元素除以相同的值 
		g[r][i]/=g[r][c]; //到过来算是因为没有存储使其变为1的约数，顺序算会导致本身变为1，则后面的元素都会除以1 
		for(int i=r+1;i<n;i++)  //在把该列该行下面的所有元素消为0
		    if(fabs(g[r][c])>eps)  //该元素不为0，则操作，否则无需操作
			for(int j=n;j>=c;j--)  //方程两边进行相同的操作 
			g[i][j]-=g[r][j]*g[i][c];  //每次操作是减掉此行第一个非0元素与上一行第一个非0元素的积，这样就能保住此行第一个非0
			                           //元素变为0 
	        r++;   //操作完一轮之后，r++，原来的最上面一行就操作完成了 
	}
	//所有行和列转换完成后，判断是否有解
	if(r<n)  //说明剩下的方程个数是小于n的，解不为一，判断是无解还是无穷多解
    {
    	for(int i=r;i<n;i++)  //r及一下的面的行，由于已经转换为阶梯行列式，所以所有系数都为0,即等式左边全为0，若等式右边
		{                      //不为0，那么左边不等于右边，则无解，否则左边等于右边，有无穷多组解
		    if(fabs(g[i][n])>eps) 
		    return 0;  //无解 
		    else 
			return 2;	//无穷多解 
		}
	        
	}
	else   //r==n  有唯一解，从下往上推出每个未知数
	{
		for(int i=n-1;i>=0;i--)  //从下往上遍历每一行,当前这一行减去下一行与当前这一行的相应系数，便于消除一个未知数 
		   for(int j=i+1;j<n+1;j++)
		   g[i][n]-=g[j][n]*g[i][j];  //等式左右两边进行相同操作
		                              //如两个方程 x1+3*x2=d1  x2=d2,那么 将1式减去2式*3则可以得出x1答案，因为是阶梯行列式
									  //所以需要依次向上推，最后一行无需推，最后一行由于先前计算阶梯矩阵已经计算出结果了 
		return 1;  //唯一解 
    }    
	      
}
int main(int argc, char** argv) {
	cin>>n;
    for(int i=0;i<n;i++)
       for(int j=0;j<n+1;j++)
       scanf("%lf",&g[i][j]);
       int t=gauss();
       if(t==1)
       {
       	 for(int i=0;i<n;i++)
       	 {
       	     printf("%.2lf\n",g[i][n]);
		 }
	   }
	   else if(t==2)
	   puts("Infinite group solutions");
	   else puts("No solution");
	return 0;
}