首页 算法 数学

# 高斯消元法

高斯消元法 是数学上线性代数中的一种算法，可以把 矩阵 转换成 行阶梯型矩阵 。

高斯消元 可以用来找到形如下列的 多元线性方程组 的解：

{a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2a31x1+a32x2+⋯+a3nxn=b3⋮am1x1+am2x2+⋯+amnxn=bm\left\{\begin{matrix} a_{11}x_1 + a_{12}x_2 + \dots +a_{1n}x_n=b_1 \\ a_{21}x_1 + a_{22}x_2 + \dots +a_{2n}x_n=b_2 \\ a_{31}x_1 + a_{32}x_2 + \dots +a_{3n}x_n=b_3 \\ \vdots \\ a_{m1}x_1 + a_{m2}x_2 + \dots +a_{mn}x_n=b_m \\ \end{matrix}\right.

对于上述方程组，我们的目标是求出所有的 未知数x 。

# 初等行列变换

我们如果将所有的系数拿出来就会形成一个系数矩阵，高斯消元法对矩阵做初等行列变换。

初等行列变换 以三种基本变换形式构成：

  1. 把某一行乘以一个非 0 的数；
  2. 交换某两行；
  3. 把某行的若干倍加到另外一行；

这三种变换都是 等价变换 。既，变换后不会对原方程的解产生任何影响

使用初等行列变换后，最理想的情况是，我们把矩阵转换成了阶梯型矩阵：

a11x1+a12x2+⋯+a1nxn=b1a22x2+⋯+a2nxn=b2⋮amnxn=bm\begin{array}{c} a_{11}x_1 + a_{12}x_2 + \dots +a_{1n}x_n & = & b_1 \\ a_{22}x_2 + \dots +a_{2n}x_n & = & b_2 \\ \vdots \\ a_{mn}x_n & = & b_m \\ \end{array}

# 如何求解

当变成如上形式之后，我们发现最后一项只有xmnx_{mn} 一个未知数，那么就可以直接求出来啦。

再看倒数第二行，只有两个未知数，既xm−1n−1x_{m-1n-1} 和mmnm_{mn} ，而其中的xmnx_{mn} 已经被求出来了，所以在倒数第二行中我们又求出了xm−1n−1x_{m-1n-1} 。以此类推，当计算完第一行之后，我们就求出了所有的 未知数x 。

从上面的求解过程中来看，显然会有三种结果：

  1. 转换后是完美阶梯形，存在唯一解；
  2. 转换后求解过程中出现了左右明显矛盾，如：左面为 0，右面不为 0。此时无解。
  3. 转换后求解过程中，某个等式出现 0=0 的情况，说明这行的等式可以取任意数，那么也就是最终方程组会有无穷多的解。

# 举个例子

当前有如下线性方程组：

{1x1+2x2−1x3=−62x1+1x2−3x3=−9−1x1−1x2+2x3=7\left\{\begin{matrix} 1x_1 + 2x_2 - 1x_3=-6 \\ 2x_1 + 1x_2 -3x_3=-9 \\ -1x_1 -1x_2 + 2x_3=7 \\ \end{matrix}\right.

转换成去掉系数的矩阵：

[12−1−621−3−9−1−127]\begin{bmatrix} 1& 2 & -1 & -6 \\ 2 & 1 & -3 \ & -9 \\ -1 & -1 & 2 & 7 \\ \end{bmatrix}

# 高斯消元法

  1. 枚举每一列。找到当前列（目前 当前列是第一列 ），绝对值最大的那一行；

  2. 将这一行换到最上面。得到：

    [21−3−912−1−6−1−127]\begin{bmatrix} 2 & 1 & -3 \ & -9 \\ 1 & 2 & -1 & -6 \\ -1 & -1 & 2 & 7 \\ \end{bmatrix}

  3. 将这行的第一个数变成 1 。也就是 2*(1/2)=1 ，那么该行所有的数都要乘上 1/2 。所以，我们就得到了：

    [10.5−1.5−4.512−1−6−1−127]\begin{bmatrix} 1& 0.5& -1.5 \ & -4.5 \\ 1& 2 & -1 & -6 \\ -1& -1 &2 & 7 \\ \end{bmatrix}

  4. 将下面所有行的当前列变成 0。那么对于第二行，它的当前列是 1，我们只要用它减去第一行，就能让它变成 0，于是，计算过后变成：

    [10.5−1.5−4.501.50.5−1.5−1−127]\begin{bmatrix} 1 & 0.5 & -1.5 & -4.5 \\ 0 & 1.5 & 0.5 & -1.5 \\ -1 & -1 & 2 & 7 \\ \end{bmatrix}

  5. 对于第三行，我们只要让它加上第一行就能把它的当前列变成 0。于是，计算后变成：

    [10.5−1.5−4.501.50.5−1.50−0.50.52.5]\begin{bmatrix} 1 & 0.5 & -1.5 & -4.5 \\ 0 & 1.5 & 0.5 & -1.5 \\ 0 & -0.5 & 0.5 & 2.5 \\ \end{bmatrix}

  6. 如此，第一轮循环就完成了，并且我们固定了第一列，剩下的方程组变成了：

    [01.50.5−1.50−0.50.52.5]\begin{bmatrix} 0 & 1.5 & 0.5 & -1.5 \\ 0 & -0.5 & 0.5 & 2.5 \\ \end{bmatrix}

  7. 继续重复步骤 1~5。目前两列中绝对值最大的是当前的第一列，将它的当前行变成 1，既， 1.5/1.5=1 ，那么该行的所有值都要 除以1.5 ，计算后得到：

    [0113−10−0.50.52.5]\begin{bmatrix} 0 & 1 & \frac{1}{3} & -1 \\ 0 & -0.5 & 0.5 & 2.5 \\ \end{bmatrix}

  8. 再把后面所有行中当前列变成 0，那么只需要让下一行加上当前行的一半： -0.5+1*1/2=0 ，那么下一行所有的数都如此做，计算后得到：

    [0113−100232]\begin{bmatrix} 0 & 1 & \frac{1}{3} & -1 \\ 0 & 0 & \frac{2}{3} & 2 \\ \end{bmatrix}

  9. 至此，已经把所有的行都完成了转换，最后得到梯形矩阵：

    [10.5−1.5−4.50113−100232]\begin{bmatrix} 1 & 0.5 & -1.5 \ & -4.5 \\ 0 & 1 & \frac{1}{3} & -1 \\ 0 & 0 & \frac{2}{3} & 2 \\ \end{bmatrix}

    对应方程组是：

    {1x1+0.5x2−1.5x3=−4.51x2+13x3=−123x3=2\left\{\begin{array}{lcl} 1x_1 + 0.5x_2 - 1.5x_3=-4.5 \\ 1x_2 + \frac{1}{3}x_3=-1 \\ \frac{2}{3}x_3=2 \\ \end{array}\right.

  10. 剩下的就是从下到上，依次求出所有的 未知数x 即可。对于最后一个方程，两面同时乘以32\frac{3}{2} ，得到x3=3x_3=3 。

    [10.5−1.5−4.50113−10013]\begin{bmatrix} 1 & 0.5 & -1.5 \ & -4.5 \\ 0 & 1 & \frac{1}{3} & -1 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix}

  11. 对于倒数第二行，只需要让它减去 最后一行乘以13\frac{1}{3} ，计算后得到：

    [10.5−1.5−4.5010−20013]\begin{bmatrix} 1 & 0.5 & -1.5 \ & -4.5 \\ 0 & 1 & 0 & -2 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix}

  12. 对于第一行，先消去第二个数 0.5 ，我们用它减去 第二行乘以 0.5 ，得到：

    [10−1.5−3.5010−20013]\begin{bmatrix} 1 & 0 & -1.5 \ & -3.5 \\ 0 & 1 & 0 & -2 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix}

  13. 对于第一行，再消去 -1.5 ，只要让它加上 第三行乘以 1.5 ，最后得到：

    [1001010−20013]\begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & -2 \\ 0 & 0 & 1 & 3 \\ \end{bmatrix}

# Code

输入格式
第一行包含整数 n。
接下来 n 行，每行包含 n+1 个实数，表示一个方程的 n 个系数以及等号右侧的常数。

输入样例：
3
1.00 2.00 -1.00 -6.00
2.00 1.00 -3.00 -9.00
-1.00 -1.00 2.00 7.00

输出样例：
1.00
-2.00
3.00

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
const int N = 110;
// 浮点数存储有误差，不能直接判 0，需要判断是否小于一个极小数。
const double eps = 1e-8;
int n;
double a[N][N]; // 存储矩阵
int gauss() {
    int c, r; //c 表示列，r 表示行
    for (c = 0, r = 0; c < n; c ++ ) {
        int t = r;
        // 找到绝对值最大的那一行。步骤 1；
        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;
        // 把这一行 换到 未处理完成行 中 的 第一行。步骤 2；
        for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);
        // 把这行的第一个数变成 1，既让该行所有的数都除以第一个数
        // 第一个数就变成了 1。这里需要注意必须从后向前算，最后一个更新第一个数。
        for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];
        // 用当前行将下面所有的列消成 0；步骤 4，5；
        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 ++ ;
    }
    if (r < n) {
        for (int i = r; i < n; i ++ ) {
            if (fabs(a[i][n]) > eps) {
                return 2; // 无解；
            }
        }
        return 1; // 无穷解；
    }
    for (int i = n - 1; i >= 0; i -- ) {
        for (int j = i + 1; j < n; j ++ ) {
            a[i][n] -= a[i][j] * a[j][n];
        }
    }
    return 0;  // 唯一解；
}
int main() {
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) {
        for (int j = 0; j < n + 1; j ++ ) {
            scanf("%lf", &a[i][j]);
        }
    }
    int t = gauss();
    if (t == 0) {
        for (int i = 0; i < n; i ++ ) {
            if (fabs(a[i][n]) < eps) a[i][n] = 0;  // 去掉输出 -0.00 的情况
            printf("%.2lf\n", a[i][n]);
        }
    }
    else if (t == 1) puts("Infinite group solutions");
    else puts("No solution");
    return 0;
}
  
 更新于   阅读次数  次
请我喝[茶]~(￣▽￣)~*
A Cat Without Sugar 微信支付
微信支付
A Cat Without Sugar 支付宝
支付宝