学习笔记 - 高斯-约旦消元法

用法

• 输入：一个矩阵，表示一个线性方程组。例如：
  $\begin{bmatrix}2&-1&1&1\\4&1&-1&5\\1&1&1&0\end{bmatrix}$，
  表示
  $\begin{cases}2x_1 - x_2 + x_3 = 1\\4x_1 + x_2 - x_3 = 5\\x_1 + x_2 + x_3 = 0\end{cases}$。

• 输出：若方程有唯一解，输出唯一的一组解。如：
  $\begin{cases}x_1 = 1\\x_2 = 0\\x_3 = -1\end{cases}$。
  算法可以判断无解/有无限组解。

• 效率：时间：$O(N^3)$，空间：$O(N^2)$。

---

做法

step1: 消元

以 $\begin{bmatrix}2&-1&1&1\\4&1&-1&5\\1&1&1&0\end{bmatrix}$ 为例，
首先将一至三行中 $x_1$ 项系数绝对值最大的一行（即第二行）移到第一行，得：

$\begin{bmatrix}4&1&-1&5\\2&-1&1&1\\1&1&1&0\end{bmatrix}$。

接着将其他方程中 $x_1$ 项的系数化为0。令第二条方程加上第一条的 $-\frac{1}{2}$ 倍，得：

$\begin{bmatrix}4&1&-1&5\\0&-\frac{3}{2}&\frac{3}{2}&-\frac{3}{2}\\1&1&1&0\end{bmatrix}$。

令第三条方程加上第一条的 $-\frac{1}{4}$ 倍，得：

$\begin{bmatrix}4&1&-1&5\\0&-\frac{3}{2}&\frac{3}{2}&-\frac{3}{2}\\0&\frac{3}{4}&\frac{5}{4}&-\frac{5}{4}\end{bmatrix}$。

接着消第二个元。

$\begin{bmatrix}4&0&0&4\\0&-\frac{3}{2}&\frac{3}{2}&-\frac{3}{2}\\0&0&\frac{17}{4}&-\frac{17}{4}\end{bmatrix}$。

消第三个元。

$\begin{bmatrix}4&0&0&4\\0&-\frac{3}{2}&0&0\\0&0&\frac{17}{4}&-\frac{17}{4}\end{bmatrix}$。

至此，消元完毕。得：

$\begin{cases}4x_1=4\\-\frac{3}{2}x_2=0\\\frac{17}{4}x_3=-\frac{17}{4}\end{cases}$。

step2: 系数化为1

系数化一，将矩阵化为单位矩阵，得：

$\begin{bmatrix}1&0&0&1\\0&1&0&0\\0&0&1&-1\end{bmatrix}$。

所以，原方程组的解就是：

$\begin{cases}x_1 = 1\\x_2 = 0\\x_3 = -1\end{cases}$。

step3: 检查无解/无穷解

如果在消元时，碰到某一元系数为0，就先跳过。到最后再来看，只要有一个出现 $bx=0(b\ne0)$ 就是无解，如果都是 $0x=0$，就是无限组解。

可以发现，对于每个元，需要 $O(N^2)$ 的时间来对每一项加减，来消元。这就是其时间复杂度由来。

---

例题

P3389 【模板】高斯消元法

模板题。不需要判断无解或无穷解，只要碰到某元系数为0，即可跳出。

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

double matrix[100][101];

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    // 读入
    int N(0);
    cin >> N;
    for (int i(0); i != N; ++i)
        for (int j(0); j <= N; ++j)
            cin >> matrix[i][j];
    // 加减消元
    for (int i(0); i != N; ++i) {
        // 找到系数绝对值最大的行
        int mx(i);
        for (int j(i + 1); j != N; ++j)
            if (fabs(matrix[j][i]) > fabs(matrix[mx][i]))
                mx = j;
        if (i != mx)
            swap(matrix[i], matrix[mx]);
        if (!matrix[i][i]) { // 无解/无限组解
            cout << "No Solution" << endl;
            return 0;
        }
        // 加减
        for (int j(0); j != N; ++j) {
            if (j != i) {
                double quot(matrix[j][i] / matrix[i][i]);
                for (int k(0); k <= N; ++k) {
                    matrix[j][k] -= matrix[i][k] * quot;
                }
            }
        }
    }
    // 系数化为1 & 输出
    for (int i(0); i != N; ++i) {
        cout << setprecision(2) << fixed << matrix[i][N] / matrix[i][i] << endl;
    }

    return 0;
}

P2455 [SDOI2006]线性方程组

其实也是模板，比上一个模板更细一些。注意细节！

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

double matrix[50][51];

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);

    int N(0);
    cin >> N;
    for (int i(0); i != N; ++i)
        for (int j(0); j <= N; ++j)
            cin >> matrix[i][j];
    int cur(0); // 这是一个细节！找系数绝对值最大的行时，范围是 [cur, N)，关注 cur 的变化
    for (int i(0); i != N; ++i) {
        int mx(cur);
        for (int j(cur + 1); j != N; ++j)
            if (fabs(matrix[j][i]) > fabs(matrix[mx][i]))
                mx = j;
        if (!matrix[mx][i]) // 无解
            continue;
        if (cur != mx)
            swap(matrix[cur], matrix[mx]);
        for (int j(0); j != N; ++j) {
            if (j != cur) {
                double quot(matrix[j][i] / matrix[cur][i]);
                for (int k(0); k <= N; ++k) {
                    matrix[j][k] -= matrix[cur][k] * quot;
                }
            }
        }
        ++cur;
    }
    if (cur != N) {
        while (cur != N)
            if (matrix[cur++][N]) { // 无解 只要有一个无解，无限组解也是无解
                cout << -1 << endl; return 0;
            }
        cout << 0 << endl; return 0; // 无限组解
    }
    for (int i(0); i != N; ++i) {
        cout << 'x' << i + 1 << '=' << setprecision(2) << fixed << matrix[i][N] / matrix[i][i] << endl;
    }

    return 0;
}