高斯消元 笔记

本文原在 2024-07-23 19:33 发布于本人洛谷博客。

一、高斯消元

如果有以下方程组，可以尝试用本文给出的方法求解：

\[\left\{\begin{matrix} 4x_1+0x_2+6x_3=28\\ 5x_1+1x_2+4x_3=22 \\ 0x_1+5x_2+1x_3=9 \end{matrix}\right. \]

把它的系数提出来，得到以下矩阵：

\[\begin{Bmatrix} 4 & 0 & 6 & | & 28\\ 5 & 1 & 4 & | & 22\\ 0 & 5 & 1 & | & 9 \end{Bmatrix} \]

基本思路是：通过更换各个方程的位置，并用加减消元法使得 \(i\) 式变为 \(\sum_{j=i}^na_{i,j}x_j=b_j\) 的形式，这样子 \(n\) 式就是 \(a_{n,n}x_n=b_n\)，用代入法逐步回带即可。

看不懂上面的式子看下面的例子，可以看出从下往上一直代入即可解出所有解：

\[\left\{\begin{matrix} x_1+x_2+x_3+x_4+x_5=5 \\ x_2+x_3+x_4+x_5=4 \\ x_3+x_4+x_5=3 \\ x_4+x_5=2 \\ x_5=1 \end{matrix}\right. \]

再说解这个方程：

• 从 \(1\) 式往下找，找到第一个有系数的 \(x_1\)，并交换到 \(1\) 式的位置。

\[\left\{\begin{matrix} 4x_1+0x_2+6x_3=28\\ 5x_1+1x_2+4x_3=22\\ 0x_1+5x_2+1x_3=9 \end{matrix}\right. \begin{Bmatrix} 4 & 0 & 6 & | & 28 \\ 5 & 1 & 4 & | & 22 \\ 0 & 5 & 1 & | & 9 \end{Bmatrix} \]

• 将 \(1\) 式乘上 \(\frac{a_{2,1}}{a_{1,1}}\) 得到 \(1'\) 式（下面未给出），即 \(\frac{5}{4}\)，在令 \(2\) 式等于 \(2\) 式减 \(1'\) 式：

\[\left\{\begin{matrix} 4x_1+0x_2+6x_3=28\\ 0x_1+1x_2-\frac{7}{2}x_3=-13\\ 0x_1+5x_2+1x_3=9 \end{matrix}\right. \begin{Bmatrix} 5 & 1 & 4 & | & 28 \\ 0 & 1 & -\frac{7}{2} & | & -13 \\ 0 & 5 & 1 & | & 9 \end{Bmatrix} \]

• 对接下来每个 \(x_i\) 和 \(i\) 式做相同操作：

\[\left\{\begin{matrix} 5x_1+1x_2+4x_3=22\\ 0x_1-1x_2+\frac{7}{2}x_3=8\\ 0x_1+0x_2+\frac{37}{2}x_3=74 \end{matrix}\right. \begin{Bmatrix} 5 & 1 & 4 & | & 22 \\ 0 & -1 & \frac{7}{2} & | & 8 \\ 0 & 0 & \frac{37}{2} & | & 74 \end{Bmatrix} \]

• 解出 \(x_3=4\)，代入 \(2\) 式，解出 \(x_2=1\)，将全部代入 \(1\) 式，解出 \(x_1=1\)。

高斯消元就是对这个矩阵进行如上操作，发现上面操作可以非常机械化的解方程，易于程序实现。

const double eps = 1e-7; // 绝对值小于 10^(-7) 的数我们就当它是 0 了
void gauss() {
    for (int i = 1; i <= n; i++) {
        int p = 0;
        for (int j = i; j <= n; j++)
            if (abs(a[j][i]) > eps) {
                p = j;
                break;
            }
        if (!p) {
            cout << "No Solution";
            return;
        }
        swap(a[i], a[p]);
        for (int j = i + 1; j <= n; j++) {
            double x = a[j][i] / a[i][i];
            for (int k = i; k <= n + 1; k++)
                a[j][k] -= a[i][k] * x;
        }
    }
    for (int i = n; i; i--) {
        a[i][i] = a[i][n + 1] / a[i][i];
        for (int j = i - 1; j; j--) {
            a[j][n + 1] -= a[j][i] * a[i][i];
            a[j][i] = 0;
        }
    }
    for (int i = 1; i <= n; i++)
        cout << fixed << setprecision(2) << a[i][i] << "\n";
}

无穷多个解：出现 \(a_{i,1}\sim a_{i,n}=0\)，\(b_i=0\) 的情况。

无解：出现 \(a_{i,1}\sim a_{i,n}=0\)，\(b_i\neq 0\) 的情况。

二、高斯-约旦消元

无需回带，直接把每一行都用加减消元消到只剩 \(a_{i,i}x_i=b_i\)。

找最大的系数来作为被减式是为了减少误差。

void gauss() {
    for (int i = 1; i <= n; i++) {
        int p = i;
        for (int j = 1; j <= n; j++) {
            if (abs(a[j][j]) > eps and j < i)
                continue;
            if (abs(a[j][i]) > abs(a[p][i]))
                p = j;
        }
        swap(a[i], a[p]);
        if (abs(a[i][i]) <= eps) continue;
        for (int j = 1; j <= n; j++) {
            if (i == j) continue;
            double x = a[j][i] / a[i][i];
            for (int k = i; k <= n + 1; k++)
                a[j][k] -= a[i][k] * x;
        }
    }
}

如果出现 \(a_{i,i}=0\)，\(b_i\neq 0\)，无解。

如果出现 \(a_{i,i}=0\)，\(b_i=0\)，无穷解。