高斯消元解线性方程组----C++实现

高斯消元$O(n^3)$

对于$n$个未知数$n$个方程构成的线性方程组：
$$\{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ \cdots \\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n\end{array}$$
其对应的增广矩阵为：
$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& b_1\\ a_{21}& a_{22}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}& b_n\end{array}\right]$$
用初等行列变化，我们要把这个方程的增广矩阵消成上三角形（唯一解），或者上阶梯型（无解或无穷解）。

矩阵有3种初等行列变化：

1. 把某一行乘一个非0的数 (方程的两边同时乘上一个非0数不改变方程的解)；
2. 交换某两行 (交换两个方程的位置)；
3. 把某行的若干倍加到另一行上去 （把一个方程的若干倍加到另一个方程上去）；

算法步骤：

从左往右枚举每一列$c$：

1. 找到当前列绝对值最大的元素所在的行；
2. 把这一行换到最上面。这里的最上面指的是还未被消成阶梯型的行中的最上面，不是整个矩阵的最上面一行；
3. 将该行乘以一个数，使得第一个元素（非阶梯型算起的第一个元素，是矩阵的第$c$列）变成1；
4. 用这个1，利用该行把该列下面的元素都消成0，注意一行的元素要一起变化。

上面循环执行完后，顺利的话，我们会得到一个上三角形：
$$\left[\begin{array}{ccccc}1& a_{12}^{\prime }& \cdots & a_{1n}^{\prime }& b_1^{\prime }\\ 0& 1& \cdots & a_{2n}^{\prime }& b_2^{\prime }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& b_n^{\prime }\end{array}\right]$$
上三角形对应的唯一解这么求：

1. 根据最后一个方程得出$x_n=b_n^{\prime }$;
2. 将$x_n$代入第$n-1$个方程，得出$x_{n-1}=b_{n-1}^{\prime }-a_{n-1,n}^{\prime }{x}_n$;
3. 不断向上回代直到第一个方程，得到$x_1$，解毕。

但有时候不能得到上三角形，只能得到一个类似阶梯型。比如我们现在的$n=4$的增广矩阵消到这个情况：
$$\left[\begin{array}{cccccc}1& 3& 2& 4& |& 5\\ 0& 1& 0& 3& |& 4\\ 0& 0& 0& 1& |& 4\\ 0& 0& 0& 4& |& 3\end{array}\right]$$
此时我们准备消第3列。但是发现，这一列剩下的元素（$(3,3),(4,3)$）中绝对值最大的元素是0，也就是说这一列全是0。那不管我们怎么搞，这一列都整不出来数字1。所以对角线位置$(3,3)$上不能是1了。没办法，只能往右边走一列。现在我们走到第4列，这一列的两个元素1,4，有非零绝对值的，所以可以正常消。消完之后得到：
$$\left[\begin{array}{cccccc}1& 3& 2& 4& |& 5\\ 0& 1& 0& 3& |& 4\\ 0& 0& 0& 1& |& 3/4\\ 0& 0& 0& 0& |& 15/4\end{array}\right]$$
记住，只要中间有一次整列都是0，那最后底下就会多出一行系数全是0的方程（比如这里在消第3列的时候出现全是0的情况，所以最后一行对应出现一个系数全为0的方程）。我们要检查底下的所有系数全是0的方程。如果

1. 出现一次诸如$0=15/4$，那方程直接无解。因为这个方程永远无法满足；
2. 如果所有方程都是$0=0$，那说明有无效不矛盾方程，方程有无穷解。

代码中，最复杂的是初等行列变化怎么实现，仔细看是不难的：

例题：高斯消元解线性方程组

#include <iostream>
#include <cmath>
using namespace std;

const int N = 110;
const double eps = 1e-6;

int n;
double a[N][N];

int gauss()
{
    int c, r; // column, row
    for (c = 0, r = 0; c < n; c ++)
    {
       int t = r;   // 从对角线元素开始往下遍历这一列
       for (int i = t; i < n; i ++) // 找到该列绝对值最大的元素所在的行号
           if (fabs(a[i][c]) > fabs(a[t][c]))
               t = i;

       if (fabs(a[t][c]) < eps) continue; // 最大绝对值是0，那么这一列剩下的元素全是0，不用管这列。这个地方导致后面不能r ++，也意味着底部将会增加一个系数全0的方程

       for (int i = c; i <= n; i ++) swap(a[t][i], a[r][i]); // 把最大绝对值元素所在的行换到未处理行的最上面(即当前要处理的的第r行)
       for (int i = n; i >= c; i --) a[r][i] /= a[r][c];  // 把现在的第r行的数字全部除以一个系数，使得左上角a[r][c]变成1
       for (int i = r + 1; i < n; i ++) // 把当前列下的所有数都消成0,要对应两行元素一起变化
           if (fabs(a[i][c]) > eps) //已经是0的就不用操作了，省点计算
               for (int j = n; j >= c; j --)
                   a[i][j] -= a[r][j] * a[i][c];

       r ++;
    }

    // 上面步骤走完之后，矩阵a[][]扣掉增广的最后一列系数以外，剩下的已经是个上三角阵或者阶梯阵
    if (r < n) // 说明有效的方程个数小于n,那要么无穷解，要么无解
    {
        for (int i = r; i < n; i ++){
            if (fabs(a[i][n]) > eps) // a[i][n] = b_i不等于0
                return 2;  // 无解
        return 1;  // 都是0 = 0的方程，无穷解
        }
    }

    // 唯一解，从下往上回代，得到方程的解
    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; j ++)
            cin >> a[i][j];

    int t = gauss();

    if (t == 0)
        for (int i = 0; i < n; i ++) printf("%.2f\n", a[i][n]);
    else if (t == 1) puts("Infinite group solutions");
    else puts("No solution");

    return 0;
}