概述

高斯-约旦消元是一种方便计算线性方程组的方法。下面我们介绍它的原理及实现。

原理

思想

各位还记得我们是怎么解多元一次方程的吗？核心思想就是消元二字。我们通过对方程进行等价变换来消元，再把得到的解回代解出所有未知数的解。这就是我们接下来要做的事情。

高斯消元

我们尝试把一个线性方程组写成矩阵的形式，比如
\begin{cases} 2x +y -z =8 \\ -3x -y +2z =-11 \\ -2x +y +2z =-3 \end{cases}
可以写成
\left[ \begin{array}{rrr|r} 2&1&-1&8 \\ -3&-1&2&-11 \\ -2&1&2&-3 \end{array} \right]
对于这个矩阵，我们尝试把左边的部分转化为一个行梯阵势。首先对于当前正在处理的行i，我们需要找到一个r>i且绝对值最大的a_{r, i}，交换r行和当前行，再继续做处理，这样可以提高数值稳定性。
我们对于正在处理的行i，使下面的每一行j的a_{j, i}变为0，即整个方程j减去a_{j, i} / a_{i, i}倍的方程i。这样，对于每一行，下面的行的第一个非0项就在第i项之后了。如果这个方程组有解，最后一个方程可以解出最后一个未知数的值，再回代求所有未知数即可。

高斯-约旦消元

上面讲述了高斯消元的过程，那么高斯-约旦消元与高斯消元有何不同呢？就是下面要介绍的内容了。
如果我们能通过类似的等价变换得到一个除主对角线以外的地方全0的矩阵，那就不用回代解未知数了，岂不美哉？。我们考虑如何用变换实现这一目标。
如果我们处理的地方不只是i行下面的每一行，而是除了i行以外的每一行。我们尝试把i列上除了a_{i, i}以外的元素全部变为0。因为在计算到该行的时候，i以前的系数全部变为0了，因此只需要对所有行的i以后的系数做变化。这样计算一番后，矩阵中留下的系数就只有a_{i, i}了，也就是说，方程组变成了a_{i, i}x_i = a_{i, n + 1}，做一下除法就求出了未知数。
这种方法的计算比高斯消元要多一些，但是省去了回代的过程，代码也短些。但是不可避免的就是复杂度仍然是O(n^3)。

代码

本代码可以通过【P3389】【模板】高斯消元法 – 洛谷一题。

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>

#include <algorithm>

const int MAXN = 105;
const double EPS = 1e-10;

int n;
double mat[MAXN][MAXN];
bool fail;

inline void gauss_jordan() {
    for(int i = 1; i <= n; i++) {
        int r = i;
        for(int j = i + 1; j <= n; j++) {
            if(std::fabs(mat[r][i]) < std::fabs(mat[j][i])) {
                r = j;
            }
        }
        if(r != i) {
            for(int j = 1; j <= n + 1; j++) {
                std::swap(mat[i][j], mat[r][j]);
            }
        }
        if(fabs(mat[i][i]) < EPS) {
            fail = true;
            return;
        }
        for(int j = 1; j <= n; j++) {
            if(j != i) {
                double t = mat[j][i] / mat[i][i];
                for(int k = i + 1; k <= n + 1; k++) {
                    mat[j][k] -= mat[i][k] * t;
                }
            }
        }
    }
    for(int i = 1; i <= n; i++) {
        mat[i][n + 1] /= mat[i][i];
    }
} 

int main() {
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n + 1; j++) {
            scanf("%lf", &mat[i][j]);
        }
    }
    gauss_jordan();
    if(fail) {
        puts("No Solution");
    } else {
        for(int i = 1; i <= n; i++) {
            printf("%.2lf\n", mat[i][n + 1]);
        }
    }
    return 0;
}

参考资料

• Gaussian elimination – Wikipedia
• 高斯消元法 – 维基百科，自由的百科全书
• 高斯-若尔当消元法 – 维基百科，自由的百科全书
• 《算法竞赛入门经典 训练指南》，刘汝佳，陈锋著，清华大学出版社，9787302291077
• [原创]选主元的高斯-约当（Gauss-Jordan）消元法解线性方程组/求逆矩阵 – 编码无悔 / Intent & Focused