高斯消元解线性方程组

基础知识

高斯消元可以通过初等行列变化把 增广矩阵 转换成 阶梯型矩阵，进而求解 n 个线性方程组的解，其时间复杂为 O(n^3)

初等行列变换：(对一个方程组进行以下三个操作不会影响方程的解)

• 把某一行乘一个非0的数（方程左右两边同乘）
• 交换某两行的位置
• 把某一行的若干倍加到另一行上

例如线性方程组为：

a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + … + a_1nx_n = b₁
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + … + a_2nx_n = b₂
…
a_n1x₁ + a_n2x₂ + a_n3x₃ + … + a_nnx_n = b_n

则增广矩阵为：

a11  a12  a13 ...  a1n  b1
a21  a22  a23 ...  a2n  b2
              .       
              .
              .
a21  a22  an3 ...  ann  bn

经过三种初等行列变换，上面的增广矩阵就可以形成一个阶梯型矩阵：

a11  a12  a13 ...  a1n  b1
     a22  a23 ...  a2n  b2
              .       
              .
              .
                    ann  bn

---

算法步骤

枚举每一列 c：
① 找到列 c 绝对值最大的一行
② 将该行换到最上面
③ 将该行第一个数变成 1 （即第 c 列的数变成 1）
④ 将下面所有行第 c 列消成 0

每枚举一列，都能将一行换到最上面(绝对值最大为 0 除外)，换到最上面说明当前这行已经固定了，所以，下次枚举某一列时，将绝对值最大的那一行换到最上面时(即②步骤)，并不是换到第一行，而是换到没有固定的行的最上面，同样，步骤①也是在未固定的行中找绝对值最大的。

答案有三种解：
如果经过以上步骤得到的是一个完美阶梯型矩阵（即得到 n 个唯一方程，即每个方程不会被另一个方程表示出来），说明有唯一解。

如果转换后出现 0==0 的方程，说明这个方程可以被其他方程表示出来，说明有 多组解。(n元方程组要想有唯一解，需要有n个唯一的方程)

如果转换后出现 0==非0 的方程，说明无解。

---

举个例子模拟一下算法

若方程组为：
x₁ + 2x₂ - x₃ = -6
2x₁ + x₂ - 3x₃ = -9
-x₁ - x₂ + 2x₃ = 7

则增广矩阵为：

1  2  -1  -6
2  1  -3  -9
-1  -1  2  7

枚举第 1 列：

① 找到第 1 列绝对值最大的一行

② 将该行换到最上面

③ 将该行第一个数变成 1（即第 1 列）
注：橙色字体为已经固定的行

④ 将下面所有行第 1 列消成 0

枚举第 2 列：
① 找到第 2 列绝对值最大的一行

② 将该行换到最上面

该行其实已经在未固定行中的最上面了

③ 将该行第一个数变成 1（即第 2 列）

④ 将下面所有行第 c 列消成 0

枚举第 3 列：
① 找到第 3 列绝对值最大的一行

② 将该行换到最上面

该行其实已经在未固定行中的最上面了

③ 将该行第一个数变成 1 （即第 c 列的数变成 1）

④ 将下面所有行第 c 列消成 0
这一行已经是最后一行了，没有下面的行了

枚举完所有列后，我们发现这个矩阵是完美阶梯矩阵，说明是有唯一解的

确定有唯一解后，从最后一行开始向上回代，第三行直接可求出 x₃ = 3，然后将 x₃ 回代到第二行，求出 x₂ = -2，再将 x₁、x₂ 回代到第一行，求出 x₁ = 1。
综上
x₁ = 1
x₂ = -2
x₃ = 3

---

例题和代码

AcWing 883. 高斯消元解线性方程组

#include<iostream>
#include<cmath>
using namespace std;
const int N = 105;
const double eps = 1e-6; //如果 x < eps，就说明 x 是等于0的 
int n;
double a[N][N];

int gauss()
{
   int c, r; // c 是列数、r 是行数
   for(c = 0, r = 0; c < n; c++) //枚举每一列 
   {
   	  int t = r;
   	  for(int i=r; i<n; i++) // 找到第 c 列的值最大的是哪一行 
   	   if(fabs(a[i][c]) > fabs(a[t][c]))
   	    t = i;
   	    
   	  if(fabs(a[t][c]) < eps) continue; //当前这一列最大值都为 0，说明这一列所有系数都为0
	   
	   for(int i=c; i<=n; i++) swap(a[t][i],a[r][i]); //把绝对值最大的这行(即第t行)换到当前最上面去(第r行)
	   
	   for(int i=n; i>=c; i--) a[r][i] /= a[r][c]; //把第t行第一个数变成1(方程所有数同除第一个数即可)
	                                               //要先处理最后一个数，因为如果先处理第一个数，第一个数先变为1了，后面所有数就会都除1而不是第一个数
												 
       for(int i =r+1; i<n; i++) //枚举从 r+1 的每一行，把第 c 列的数全部消成 0 
	      if(fabs(a[i][c]) > eps) //如果已经 ==0 ，就不用消了
		    for(int j=n; j>=c; j--) //枚举这一行的每个数 
			  a[i][j] -= a[r][j] * a[i][c]; //每个数减去a[i][c] * a[r][j] 即可 
		  												   
		r++; //该下一行了										   	   
   }
   
   if(r < n) //说明剩下的方程数小于 n 个，需要判断是 无解 还是 有多组解 
   { //现在已经是阶梯型，所有系数都应该为0，如果等式右边a[i][n] != 0，说明无解 
   	  for(int i=r; i<n; i++)
   	    if(fabs(a[i][n]) > eps)
   	      return 2;  //无解
   	  return 1; //多组解，说明 r~n-1的这些方程是可以被 0~n 的某些方程表示出来的
		   
   }
   
   //有唯一解，从下向上回代，依次求解 
   for(int i=n-1; i>=0; i--) //从下向上枚举每一行 
     for(int j=i+1; j<n; j++) 
       a[i][n] -= a[i][j] * a[j][n];
   
   return 0;   
} 

int main()
{
   	cin >> n;
   	for(int i=0; i<n; i++)
   	 for(int j=0; j<n+1; j++)
   	   cin >> a[i][j];
   	   
   	int t = gauss();
	if(t == 0) //有唯一解 
	{
		for(int i=0; i<n; i++) printf("%.2lf\n",a[i][n]);
	}      
	else if(t == 1) printf("Infinite group solutions\n"); //有多组解
	else printf("No solution\n"); //无解 
} 

---