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

题目传送门

一、高斯消元 \(O(n^3)\)

• 通过初等行变换把增广矩阵化为阶梯型矩阵并回代得到方程的解。

• 适用于求解 包含\(n\) 个方程，\(n\) 个未知数的多元线性方程组。

例如该方程组

$ \left\{ \begin{array}{lc} a_{11}x_1 +a_{12}x_2+a_{13}x_3+......+a_{1n}x_n =b_1 \\ a_{21}x_1 +a_{22}x_2+a_{23}x_3+......+a_{2n}x_n =b_2 \\ a_{31}x_1 +a_{32}x_2+a_{33}x_3+......+a_{3n}x_n =b_3 \\ a_{41}x_1 +a_{42}x_2+a_{43}x_3+......+a_{4n}x_n =b_4 \\ a_{51}x_1 +a_{52}x_2+a_{53}x_3+......+a_{5n}x_n =b_5 \\ ... \\ a_{n1}x_1 +a_{n2}x_2+a_{n3}x_3+......+a_{nn}x_n =b_n \\ \end{array} \right. $

增广矩阵为：

\[\begin{gathered} \begin{pmatrix} a_{11} & a_{12} & ... &a_{1n} & b_1 \\ a_{21} & a_{22} & ... &a_{2n} & b_2 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n1} & a_{n2} & ... &a_{nn} & b_n \end{pmatrix} \quad \end{gathered} \]

• 接下来的所有操作都用该增广矩阵，代替原方程组

二、前置知识：初等行（列）变换

1. 方程两边同时乘上一个非\(0\)数,不改变方程的解。

2. 交换两个方程的位置。

3. 把某行的若干倍加到另一行上去,起到消元的效果。

接下来，运用初等行变换，把增广矩阵变为阶梯型矩阵。

阶梯型矩阵:

\[\begin{gathered} \begin{pmatrix} a_{11} & a_{12} & ... &a_{1n} & b_1 \\ & a_{2i} & a_{2(i+1)} &a_{2n} & b_2 \\ & & & ⋮ & ⋮ \\ & & &a_{nn} & b_n \end{pmatrix} \quad \end{gathered} \]

最后再把阶梯型矩阵从下到上回代到第一层即可得到方程的解。

三、算法步骤

. 枚举每一列\(c\)，找到当前列绝对值最大的一行

. 用初等行变换(\(2\)) 把这一行换到最上面（未确定阶梯型的行，并不是第一行）

. 用初等行变换(\(1\)) 将该行的第一个数变成 \(1\) （其余所有的数字依次跟着变化）

. 用初等行变换(\(3\)) 将下面所有行的当且列的值变成 \(0\)

实例演示：

四、完整代码

#include <bits/stdc++.h>

using namespace std;

const int N = 110;
const double eps = 1e-6; //浮点数精度控制

int n;
double a[N][N]; //系数+结果矩阵 n*(n+1)
int c, r;               //c表示是枚举的列,r表示是枚举的行

// 高斯消元求多元一次方程组
int gauss() {
    /**
     如果有唯一解，从最后一个式子开始向上逐渐消元直至得到最简上三角行列式
     */
    for (c = 0; c < n; c++) {
        // 第一步：找到c列绝对值最大的行数
        int t = r; //暂存的最大行,猴子选大王
        for (int i = r; i < n; i++)
            if (fabs(a[i][c]) > fabs(a[t][c])) t = i;

        if (fabs(a[t][c]) < eps) continue; //如果没有找到系数大于0的行，无需调整

        // 第二步：将t行换到最上方还未确定的一行
        for (int i = c; i < n + 1; i++) swap(a[t][i], a[r][i]);

        // 第三步：用初等行变换(1)将该行的第一个数变成1（其余所有的数字依次跟着变化）
        //方程两边同时除以第一个数a[r][c]，需要倒着算，不然第一个数先变1，系数就被篡改
        //后面的数字不知道该除谁了。
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c];

        //(4)用初等行变换(3)将下面所有行的当前第C列的值清为0
        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++;
    }
    /**
    * 如果是唯一解，那么每一行确定一个一个解
    * 假设有3个方程，3个未知量，那么r=0时处理第一个未知量，=1时第二个，=2时第三个，
    * =3时结束了，也就是r最终停止的位置是不确定解的位置
    * 所以如果是无解或者无穷解，r-1是最后一个确定方程解得位置，r是第一个全零行，
    * 所以判断是无解还是无穷解是从r开始的，而非r+1
    */
    // 第五步：判断是否有解并处理可消去可消系数
    if (r < n) {
        for (int i = r; i < n; i++)
            if (fabs(a[i][n]) > eps)
                return 2;//无解
        return 1;//无穷多组解
    }
    //第六步：用i行下面每一行来把i行对应的系数消为0，常数列对应改变
    for (int i = n - 1; i >= 0; i--)        //i行
        for (int j = i + 1; j < n; j++)     //j列，从i+1开始。至n-1，不包含最后列的方程等号右边的数字
            a[i][n] -= a[j][n] * a[i][j];   //每消一行，其实就剩下an*xn=b,其它的都是0了。

    //有唯一解
    return 0;
}

int main() {
    //优化输入
    ios::sync_with_stdio(false);
    cin >> n;
    for (int i = 0; i < n; i++)     //n*(n+1) 矩阵
        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;
}