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

前置

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

模板题

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

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

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

补充

• 它和普通的高斯消元法有什么区别吗？
  高斯消元法每次消元，最后会得到系数矩阵的上三角矩阵
  然后通过回带来得到所有解，稍微麻烦些。
  举个例子，比如得到下面的样子：
  $$A=\left[\begin{array}{ccccc}2& 4& 12& |& 10\\ 0& 8& 2& |& 5\\ 0& 0& 4& |& 20\end{array}\right]$$
  然后易得 $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;											/// 有唯一解
}