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;
}