【算法讲13：高斯约旦消元法】改进：判断无解与无穷解

前置在线性代数的课程中，我们就学过基本的矩阵及其行初等变化。根据这些初等变化，我们的老师就高速我们怎么样去进行消元，然后求解线性方程组。模板题【模板】高斯消元法 | 洛谷 P3389挺基础的蓝题，虽然是个模板但是思路很简单可以手敲。⌈\lceil⌈高斯约旦消元法⌋\rfloor⌋首先，我们获得一个系数矩阵 An×nA_{n\times n}An×n 和一个列向量 Bn×1B_{n\times 1

溢流眼泪

2851人浏览 · 2021-02-13 20:01:43

溢流眼泪 · 2021-02-13 20:01:43 发布

前置

在线性代数的课程中，我们就学过基本的矩阵及其行初等变化。
根据这些初等变化，我们的老师就高速我们怎么样去进行消元，然后求解线性方程组。

模板题

【模板】高斯消元法 | 洛谷 P3389
挺基础的蓝题，虽然是个模板但是思路很简单可以手敲。

$\lceil$ 高斯约旦消元法 $\rfloor$

首先，我们获得一个系数矩阵 $A_{n\times n}$ 和一个列向量 $B_{n\times 1}$
我们把它合并成一个增广矩阵 $A_{n\times (n+1)}$
那么这个矩阵代表什么意思呢？比如我们有一下方程组：
$\begin{cases} 1x+8y+2z=5\\ 2x+4y+12z=10\\ 1x+0y+4z=20\\ \end{cases}$
我们的增广矩阵就如下所示
$\begin{bmatrix} 1&8&2&|&5\\ 2&4&12&|&10\\ 1&0&4&|&20\\ \end{bmatrix}$
【做法】：
（1）我们按列进行处理。假设现在处理到第 $i$ 列。
为了使得答案精度更高，我们选择 该列中行号 $\ge i$ 的行中最大的值的那行 ，然后把该行和第 $i$ 行进行交换。
（2）如果某一列的所有值都为 $0$ ，那么说明该行变量是个无关变量，任意取值都满足。
（3）根据行初等变化，我们把除了第 $i$ 行的每一行都减去一定的比值，使得 该行第 $i$ 列为 $0$
（4）这样操作完，系数矩阵只有主对角线位置的元素非零。容易算出每一个变量的取值了。
【举例】：
$\begin{bmatrix} 1&8&2&|&5\\ 2&4&12&|&10\\ 1&0&4&|&20\\ \end{bmatrix}$
（1）我们找第一列最大的元素，为 $2$ ，该行和第一行交换
$\begin{bmatrix} 2&4&12&|&10\\ 1&8&2&|&5\\ 1&0&4&|&20\\ \end{bmatrix}$
（2）第二行减去第一行的 $\frac{1}{2}$ 倍，第三行减去第一行的 $\frac{1}{2}$ 倍，得到新的矩阵：
$\begin{bmatrix} 2&4&12&|&10\\ 0&6&-4&|&0\\ 0&-2&-2&|&15\\ \end{bmatrix}$
（3）找第二列的行号为 $2\sim 3$ 的最大值。由于 $6$ 最大且本身就是第二列，因此无需交换。
（4）第一行减去第二行的 $\frac{4}{6}$ 倍，第三行减去第二行的 $-\frac{2}{6}$ 倍，得到：
$\begin{bmatrix} 2&0&\frac{44}{3}&|&10\\ 0&6&-4&|&0\\ 0&0&-\frac{10}{3}&|&15\\ \end{bmatrix}$
（5）找第三列的行号为 $3\sim 3$ 的最大值。最大值为 $-\frac{10}{3}$ ，无需交换。
（6）第一行减去第三行的 $-\frac{44}{3}\cdot\frac{3}{10}$ 倍，第二行减去第三行的 $4\cdot\frac{3}{10}$ 倍。得到：
$\begin{bmatrix} 2&0&0&|&76\\ 0&6&0&|&-18\\ 0&0&-\frac{10}{3}&|&15\\ \end{bmatrix}$
我们容易得到，变量的值为：
$\begin{cases} x=38\\ y=-3\\ z=-4.5\\ \end{cases}$
完美！~~(手算错两遍咳咳咳)~~

补充

它和普通的高斯消元法有什么区别吗？
高斯消元法每次消元，最后会得到系数矩阵的上三角矩阵
然后通过回带来得到所有解，稍微麻烦些。
举个例子，比如得到下面的样子：
$\begin{bmatrix} 2&4&12&|&10\\ 0&8&2&|&5\\ 0&0&4&|&20\\ \end{bmatrix}$
然后易得 $z$ 的值，然后带回第二行得到 $y$ 的值，再带回第一行得到 $x$ 的值。
我个人感觉高斯约旦消元法更加方便。

代码

时间复杂度： $O(N^3)$

/*
 _            __   __          _          _
| |           \ \ / /         | |        (_)
| |__  _   _   \ V /__ _ _ __ | |     ___ _
| '_ \| | | |   \ // _` | '_ \| |    / _ \ |
| |_) | |_| |   | | (_| | | | | |___|  __/ |
|_.__/ \__, |   \_/\__,_|_| |_\_____/\___|_|
        __/ |
       |___/
*/
const int MAX = 1e2+50;
const ll  MOD = 1e9+7;

double aa[MAX][MAX];

bool GJ(int n){
    for(int i = 1;i <= n;++i){
        int now = i;
        double mx  = aa[i][i];
        for(int j = i + 1;j <= n;++j){
            if(aa[j][i] > mx){
                now = j;
                mx = aa[i][j];
            }
        }
        if(fabs(mx) < EPS)return false;
        if(now != i)
            for(int j = 1;j <= n + 1;++j)
                swap(aa[now][j],aa[i][j]);
        for(int j = 1;j <= n;++j){
            if(j == i)continue;
            for(int k = n + 1;k >= i;--k){
                aa[j][k] -= aa[j][i] / aa[i][i] * aa[i][k];
            }
        }
    }
    return true;
}
int main()
{
    int n;cin >> n;
    for(int i = 1;i <= n;++i)
    for(int j = 1;j <= n + 1;++j)
        scanf("%lf",&aa[i][j]);
    bool can = GJ(n);
    if(!can)puts("No Solution");
    else{
        for(int i = 1;i <= n;++i)
            printf("%.2f\n",aa[i][n+1] / aa[i][i]);
    }
    return 0;
}

更新

突然出了这么一题线性方程组 | 洛谷 P2455
除了输出唯一解外，还需要判定无解与无穷解的情况。
虽然高斯约旦消元法不大擅长判断这种无解还是无穷解的情况，但是我们仍然可以用高斯约旦消元法去做。
主要思路就是我们设一个当前行 $r$ ，如果你找不到当前列的主元(就是全是0)，那么直接跳到下一列(但是当前行不变)；否则，就正常的消元，然后当前行递增。
当你消元完毕后，如果当前行 $r < n$ ，那么表示最后几行的系数矩阵全为 $0$ 。我们需要判断这几行的列向量是否全 $0$ ，以此判断是无解还是无穷解。

改进代码

double aa[MAX][MAX];	/// 增广矩阵
double ans[MAX];
int GJ(int n){
    int r = 1;
    for(int c = 1;c <= n;++c){
        int now = r;
        double mx  = fabs(aa[r][c]);			/// 找最大系数改为找绝对值最大系数
        for(int i = r + 1;i <= n;++i){
            if(fabs(aa[i][c]) > mx){
                now = i;
                mx = fabs(aa[i][c]);
            }
        }
        if(mx < EPS)continue;					/// 跳过当前列
        for(int j = 1;j <= n + 1;++j)
            swap(aa[r][j],aa[now][j]);
        for(int i = 1;i <= n;++i){
            if(i == r)continue;
            for(int j = c + 1;j <= n + 1;++j)
                aa[i][j] -= aa[i][c] / aa[r][c] * aa[r][j];
        }
        r++;									/// 当前行变为下一行
    }
    r--;
    if(r < n){
        for(int i = r + 1;i <= n;++i)
            if(fabs(aa[i][n + 1]) > EPS)return -1;		/// 无解
        return 0;										/// 无穷解
    }
    for(int i = 1;i <= n;++i){
        ans[i] = aa[i][n + 1] / aa[i][i];
        if(fabs(ans[i]) < EPS)ans[i] = 0;				/// 注意把输出 -0改成 0
    }
    return 1;											/// 有唯一解
}