高斯消元

abensyl · 2023-12-20 21:42:35 · 题解

P3389【模板】高斯消元法。

2023.12.21 更新：更正了小错误，添加了更详细的解释。

前置线性代数知识：线性相关、初等变换、矩阵的秩、阶梯矩阵。

思路

这道题目是一个很好的考察队线性代数掌握情况的模板题，当然，你用小学数学里学的那个消元法的实质也就是线性方程组的初等变换，而线性方程组的初等变换和矩阵的初等变换在本质上是一样的，所以我觉得用矩阵的初等变换来解释更加直观，因为矩阵的秩可以直接说明方程组解的情况，并且阶梯形矩阵的概念更加直观易懂。

解线性方程组需要消元，消元的操作可以视作对于一个增广矩阵求阶梯矩阵的过程，对于下面这个方程，我们从矩阵的角度来剖析他：

\alpha_{1,1} x_1+\alpha_{1,2} x_2+\dots\alpha_{1,n} x_n=\beta_1\\ \alpha_{2,1} x_1+\alpha_{2,2} x_2+\dots\alpha_{2,n} x_n=\beta_1\\ \dots\\ \alpha_{n,1} x_1+\alpha_{n,2} x_2+\dots\alpha_{n,n} x_n=\beta_n\\ \end{cases}

（其中 x_1,x_2,\dots,x_n 为未知量，\alpha_{i,j} 为常数）

增广矩阵就可以理解为两个矩阵的组合，而对于一个 n 元 n 个方程的线性方程组，它的增广矩阵的左侧的矩阵规模为 n\times n，也就是 n 个列向量，其中第 i 个列向量的第 j 行就是 x_i 在第 j 个方程里的系数，而右侧的矩阵是一个规模为 n\times 1 的矩阵，也可以理解为一个列向量，第 i 行的数就是第 i 个方程组。

把这个矩阵用刚才线性方程组里的字母写出来，就是：

\def\arraystretch{1.5} \begin{array}{c|c} {\begin{matrix} \alpha_{1,1}&\alpha_{1,2}&\dots&\alpha_{1,n}\\ \alpha_{2,1}&\alpha_{2,2}&\dots&\alpha_{2,n}\\ \dots&\dots&\dots&\dots\\ \alpha_{n,1}&\alpha_{n,2}&\dots&\alpha_{n,n} \end{matrix}}& {\begin{matrix} \beta_1\\\beta_2\\\dots\\\beta_n\end{matrix}} \end{array} \right)

我们要做的第一步实际就是构建这个增广矩阵，构建方法非常简单，直接读入就可以了，因为题目的输入就是按照矩阵的格式读入的，我们在这里假设矩阵为 R，简记竖线左侧的矩阵为 A，简记竖线右侧的矩阵为 B，那么显然 R=\left(\begin{array}{c|c}A&B\end{array}\right)，这里，我们习惯性的称 A 为系数矩阵，R 为增广矩阵。

那么，我们知道，矩阵是可以进行初等变换的，在这里，用采用初等行变换进行消元，初等变换有以下三种：

用一非零的数 \lambda 乘以某一整行；
把一行的 \lambda 倍加到另一行；
互换任意两行的位置。

（下如需使用，按顺序简称第一种、第二种、第三种）

我们知道，如果我们可以把一个矩阵化成阶梯矩阵的形式，那么就可以直接得到方程的解，所以我们的工作就是通过初等变换，化已知矩阵为阶梯矩阵即可。

我们知道，矩阵的秩就是矩阵中线性无关的行数的极大值（列数也一样，不过显然方程组是横着写的），简记为 \operatorname{rank}{R}，那么判断有无解的方式很显然，就是判断 A、R 两个矩阵的秩以及 n 的关系就可以了，如果 \operatorname{rank}{R}\neq \operatorname{rank}{A} 或 \operatorname{rank}{R}>n 说明无解（后者在这题里不可能出现，因为行数只有 n），否则如果 \operatorname{rank}{R}<n 说明有多个解，否则有唯一解，在这里，矩阵的秩其实就可以理解为化为阶梯矩阵后不为零的列数。

可总结为：

如果系数矩阵的秩恰等于增广矩阵的秩，且等于未知数数量就是唯一解；
如果这个系数矩阵的秩恰不等于增广矩阵的秩就无解，或者增广矩阵的秩大于未知数数量也是无解；
剩下的情况就是无数解。

那么，高斯消元就是用已经算好的行去消剩下的行，具体的消除方法如下：

如果有解，则一定秩满（\operatorname{rank}{R}=n），所以一定满足主对角线上所有元素非零，主对角线以左下所有元素为零，我们从第 1 行向第 n 行逐步消元，保证 \alpha_{i,1}=\dots=\alpha_{i,i-1}=0 且 \alpha_{i,i}\neq 0，所以我们需要检查第 i 行第 i 列，如果为零，就需要考虑跟其他第 i 列不为 0 的没被处理过的行第三种初等变换，如果所有行的第 i 列均为 0，就说明秩不满，直接输出 No Solution 即可，实质就是扫一遍没处理过的行随便找一个第 i 列不为 0 的放到第 i 行去；
用第 i 行将第 i+1\sim n 行的第 i 个元素全部消掉（在进行处理第 i 行的时候，前 i-1 行肯定都已经消掉了），那么这里需要用到第二种初等变换，具体的操作中，我们需要确定 \lambda 的值，很显然，如果我们要将第 j（j>i）行的第 i 个元素消掉，那么就满足 \alpha_{j,i}+\lambda\times\alpha_{i,i}=0，显然 \lambda=-\frac{\alpha_{j,i}}{\alpha_{i,i}}；

划为有唯一解的阶梯矩阵以后，矩阵大致形态如下：

\def\arraystretch{1.5} \begin{array}{c|c} {\begin{matrix} \theta_{1,1}&\theta_{1,2}&\dots&\theta_{1,n-1}&\theta_{1,n}\\ 0&\theta_{2,2}&\dots&\theta_{2,n-1}&\theta_{2,n}\\ \dots&\dots&\dots&\dots&\dots\\ 0&0&\dots&\theta_{n-1,n-1}&\theta_{n-1,n}\\ 0&0&\dots&0&\theta_{n,n} \end{matrix}}& {\begin{matrix} \epsilon_1\\\epsilon_2\\\dots\\\epsilon_n\\\epsilon_n\end{matrix}} \end{array} \right)

有了这样一个矩阵，那么后续就容易了，首先第 n 行是一个一元一次方程组，非常好解，然后最后一行求出 x_n=\frac{\epsilon_n}{\theta_{n,n}}，将其作为常数带入倒数第二行，得到 x_{n-1}=\frac{\epsilon_n-x_n\times \theta_{n-1,n}}{\theta_{n,n}}，以此类推，可以求出整个方程组的解。

写代码时的一些注意点：

本题在消元时使用第二种初等变换，不能保证中间结果和最终答案均为整数，需要使用浮点数；
注意判断浮点数是否为零的时候需要小心精度差，判 a[i][i]=0 的时候请用 fabs(a[i][i])<eps 而不是 a[i][i]==0。否则开了 O2 优化后会因为 FMA 挂掉（详见这篇帖子）；
如果第 i 行第 i 列为零，一定要和下面的其他行进行第二种初等变换，否则程序可能输出 nan；
注意 No Solution 的拼写方式。

代码里有简要的注释，可供参考。

代码

#include <bits/stdc++.h>
const long double eps=1e-7;
using namespace std;
long double mat[105][105],res[105];
int n;
void trans(int i,int j) {
    // 初等变换 3：交换第 i 行和第 j 行
    if(i==j) return;
    swap(mat[i],mat[j]);
}
void trans(int i,int j,double k) {
    if(!i||!j) return;  
    // 初等变换 2：将第 i 行的 k 倍与第 j 行对应相加
    for(int l=1;l<=n+1;++l)
        mat[j][l]+=k*mat[i][l];
}
signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr),
    cout.tie(nullptr);
    cin>>n;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n+1;++j) 
            cin>>mat[i][j];
    for(int i=1;i<=n;++i) {
        bool flag=0;
        // 先找到一个第 i 个数字不为 0 的行将其初等变换到第 i 行
        for(int j=i;j<=n;++j)
            trans(i-1,j,-mat[j][i-1]/mat[i-1][i-1]);
        for(int j=i;j<=n;++j)
            if(!fabs(mat[j][i])<eps) {
                trans(i,j);
                flag=1;
                break;
            }
        if(!flag) {
            // 不满秩，找不到就一定没有唯一解
            cout<<"No Solution\n";
            return 0;
        }
        if(!mat[i][i]) {
            // 初等变换把 (i,i) 也变成 0 了
            cout<<"No Solution\n";
            return 0;
        }
        // 通过初等变换 2 进行消元
    } // 变为阶梯形矩阵后，从下往上求方程的解
    for(int i=n;i>=1;--i) {
        res[i]=mat[i][n+1];
        for(int j=n;j>i;--j)
            res[i]-=mat[i][j]*res[j];
        res[i]/=mat[i][i];
        // 用之前算出的结果从下向上代入方程求解
    }
    for(int i=1;i<=n;++i)
        cout<<fixed<<setprecision(2)<<res[i]<<'\n';
    // 输出答案
    return 0;
}