高斯消元

洛谷模板链接：高斯消元法

题目描述

给定一个线性方程组，对其求解

输入输出格式

输入格式： 

第一行，一个正整数 nn

第二至 n+1n+1 行，每行 n+1n+1 个整数，为a_1, a_2 cdots a_na1​,a2​⋯an​ 和 bb ，代表一组方程。

输出格式：

共n行，每行一个数，第 ii 行为 x_ixi​ （保留2位小数）

如果不存在唯一解，在第一行输出"No Solution". 

输入输出样例

输入样例#1： 

3
1 3 4 5
1 4 7 3
9 3 2 2

输出样例#1： 

-0.97
5.18
-2.39

高斯消元法是线性代数规划中的一个算法，主要用来为线性方程组求解，也就是求解多元一次方程的解。也可以求解矩阵的秩，矩阵的逆。

下面来讲一下高斯消元的过程。

首先我们回想一下解方程组的方法，我们学过用换元法，消元法来减少未知数的个数，高斯消元就是用的消元法，我们以这样一组数据为例:

1x - 2y + 3z = 6　①
4x - 5y + 6z = 12　②
7x - 8y + 10z = 21　③

但是我们在求解方程组的时候没有必要知道未知数的名字，答案只与未知数的系数有关，所以可以直接将系数存到数组中，那么方程组就变成了这个样子：

1 -2 3 6　①
4 -5 6 12　②
7 -8 10 21　③

（前

我们一项一项的去掉未知数的系数，最后我们期待得到的样子是这样的：

1 0 0 a (1x + 0y + 0z == a)
0 1 0 b (0x + 1y + 0z == b)
0 0 1 c (0x + 0y + 1z == c)

这样的矩阵对角线的值都为1（∑f[i][i] = 1,0 <= i <= n-1，i用0~n-1作为未知数系数，所有的f[i][n]是每个方程常数项的系数）,第i行第n个表示的也就是每个未知数的值．所以我们考虑如何将原矩阵转换成这个样子.

我们用i来表示正在消去第i行的系数。为了方便消去除了第i行以外其他行的系数,我们可以先把f[i][i]的系数消为1(因为最终要将矩阵消去系数直到只有对角线有值且值为1)．

可以一开始先将第i列系数最大的值交换到第一行，每次处理第i个系数最大的那个值，这样可以减小误差．

7 -8 10 21　①
4 -5 6 12　②
1 -2 3 6　③

然后将第i个的系数转化为1.

1 -8/7 10/7 3　①
4 -5 6 12　②
1 -2 3 6　③

接下来就可以通过②-①*4和③-①消去其它行的第i个的系数（此时i == 0）．

1 -8/7 10/7 3　①
0 -3/7 2/7 0　②
0 -6/7 11/7 3　③

然后接下来每一行都是这个套路．

1 -8/7 10/7 3 ①
0 1 -2/3 0　②
0 -6/7 11/7 3　③

通过①-②*(-8/7)和③-②*(-6/7)消去第一行和第三行的第i个系数（此时i == 1）．

1 0 2/3 3 ①
0 1 -2/3 0　②
0 0 1 3 ③

最后一行不需要继续化简，直接①-③*(2/3),②-③*(-2/3)来得到最后的答案．

1 0 0 1
0 1 0 2
0 0 1 3

此时，一次完整的消元就做完了．

那么我们来总结一下高斯消元的主要过程(此时正在处理第i行,0 <= i <= n-1)：

1. 在 j = i~n-1 行中选择 f[j][i] 最大的那一行与第i行交换.
2. 将第i行第i个系数消为1,后面所有的值都除以f[i][i].
3. 将除第i行以外所有行j的所有系数f[j][k]都减去f[i][k]*f[l][i],使得除了第i行,每一行第i列的系数都被消为0.
4. 重复以上步骤直到消完最后一行.

其实方程组也存在无解和无数组解的情况,为什么会无解呢?如果给出来一组这样的数据,我们想想会怎么样:

2 1 3
2 1 3

经过消元后:

1 1/2 3/2
0 0 0

这样也就是在消元过程中处理第2行时第2行的第2位是0,说明不存在方法能够确定第2个未知数的值,说明有无数组解.其实这个数据给出的两个方程是一样的,那么我们就得到了判断无数组解的方法: 判断 f[i][i]和f[i][n] 是否等于 0.若都等于0,则说明存在无数组解.

我们再把这个数据改一下,变成这个样子:

2 1 3 
2 1 2

消元后:

1 1/2 3/2
0 0 1

这样的情况和刚才那个有点像,但是这个的常数项不等于0,此时无法解出方程组.所以判断是否无解的方法也就出来了: 判断一行是否除了常数项所有系数都为0.如果有这样的一行,那么无解.

下面看一下代码注释(无解和无数组解都输出No Solution):

#include<bits/stdc++.h>
using namespace std;
const int N=100+5;
const double eps=1e-8;

int n;
double f[N][N];

int main(){
  cin >> n;
  for(int i=0;i<n;i++)
    for(int j=0;j<=n;j++)
      cin >> f[i][j];
  for(int i=0;i<n;i++){
    int pos = i;
    for(int l=i;l<n;l++)
      if(fabs(f[pos][i]-f[i][i]) <= eps) pos = l;//步骤1
    for(int k=0;k<=n;k++)
      swap(f[pos][k] , f[i][k]);//这个循环相当于是将第i行和第pos行交换了
    if(fabs(f[i][i]) <= eps){
      printf("No Solution
");//这个是判断无解的情况
      return 0;
    }
    for(int k=i+1;k<=n;k++) f[i][k] /= f[i][i];//步骤2
    for(int l=0;l<n;l++)
      if(i != l)
    for(int k=i+1;k<=n;k++)
      f[l][k] -= f[i][k]*f[l][i];//步骤3
  }
  for(int i=0;i<n;i++)
    printf("%.2lf
",f[i][n]);
  return 0;
}