[数学+模板] 高斯消元算法模板(模板)

1. 高斯消元算法+模板

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

重点： 高斯消元

高斯消元能够在 $O(n^3)$ 时间内求解 n 元一次线性方程组。也称高斯列主消元法。

具体算法步骤：

• 枚举每一列 c，循环做如下操作
  • 找到绝对值最大的一行。即，找到列主元
  • 将该行整体交换到最上面。准备进行列主消元
  • 将该行第一个数变成 1。做初等行列变换
  • 将下面所有行的第 c 列变成 0。做初等行列变换
• 处理完全部列，此时的增广矩阵即为一个上三角矩阵。最后一行的解已知，反推即可得到所有的解。

讨论解的情况：

n 元一次线性方程组解就三种情况：

• 完美阶梯型—唯一解，$a{x}_i=p$ 则 $x_i=\frac{p}{a}$ 则解 $x_i$ 唯一确定。
• 出现 0 = 非 0—无解，即有 $0x_i=p$ 则不存在这样的 $x_i$ 使等式成立，则无解
• 出现 0 = 0—无穷多解，即有 $0x_i=0$ 则存在任意多的 $x_i$ 使得等式成立，则无穷多解

高斯列主消元是数值计算中的最基础算法之一，代码简单，算法好理解。对于解的情况也可以从矩阵的秩的角度去考虑都是可以的。主要还是对二维数组的操作，是一个不错的练手项。

代码注意点：

• 浮点数绝对值函数 fabs()
• 浮点数比较，精度问题
• 在将列主元变成 1 的过程中，需要倒着将本行中的所有元素除上列主元。要是正着除，则列主元直接变成 1，就不能再继续使用列主元了，需要额外再保存一下列主元数值，比较麻烦。
• 无解的判断、列主元为 0 时的判断

模板代码：

#include <iostream>
#include <cmath>

using namespace std;

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

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

int gauss() {
    int c, r;
    for (c = 0, r = 0; c < n; ++c) {        // c列r行，遍历列
        int t = r;
        
        for (int i = r; i < n; ++i)         // 寻找列主元，拿t记录
            if (fabs(a[i][c]) > fabs(a[t][c])) 
                t = i;              
        
        if (fabs(a[t][c]) < eps) continue;  // 如果列主元为0，不必考虑，当前列全为0
        
        // 交换列主元t行与当前r行的数
        for (int i = c; i < n + 1; ++i) swap(a[t][i], a[r][i]); 
        // 当前列主元已经被交换到了r行，需要从后向前进行处理，避免主对角线元素变成1
        for (int i = n; i >= c; --i) a[r][i] /= a[r][c]; 
        
        // 列主消元
        for (int i = r + 1; i < n; ++i) 
            if (fabs(a[i][c]) > eps)
                for (int j = n; j >= c; --j) 
                    a[i][j] -= a[r][j] * a[i][c];
                    
        ++r;
    }
    
    if (r < n) {
        for (int i = r; i < n; ++i) 
            if (fabs(a[i][n]) > eps) return 2;  // 0x=1 则无解
        
        return 1;   // 0x=0 无穷多解
    }
    
    // 上三角阶梯型矩阵
    // 直接求解即可，最后一列放置结果
    for (int i = n - 1; i >= 0; --i)    
        for (int j = i + 1; j < n; ++j) 
            a[i][n] -= a[j][n] * a[i][j];
   
    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) puts("Infinite group solutions");
    else puts("No solution");
    
    return 0;
}