高斯消元（洛谷P3389）

简介

数学上，高斯消元法，是线性代数规划中的一个算法，可用来为线性方程组求解。常用于加减消元法，求出矩阵的秩，以及求出可逆方阵的逆矩阵。

实现

其实高斯消元就是模拟手算方程组的过程。

手算方程组最常用的两种方法：加减消元和代入消元。不断消元直到解出一个元，再不断回代解出其他元。

高斯消元法则是通过不断加减消元解出一个元再回代回去。

举个例子

比如下面这个方程组：

$$\begin{equation} \left\{ \begin{aligned} x+y+z&=6 \quad ① \\ 2x+3y-z&=5 \quad ②\\ x-2y+2z&=3 \quad ③ \\ \end{aligned} \right. \end{equation}$$
我们用矩阵来表示：
$$\left[ \begin{matrix} x & y & z & w \\ 1 & 1 & 1 & 6 \\ 2 & 3 & -1 & 5 \\ 1 & -2 & 1 & 0 \\ \end{matrix} \right]$$
现在我们消去 $x$：

① $\times2-$ ②得 $\ -y+3z=7\quad ②$

① $-$ ③得 $3y-z=3\quad③$

即

$$\left[ \begin{matrix} x & y & z & w \\ 1 & 1 & 1 & 6 \\ 0 & -1 & 3 & 7 \\ 0 & 3 & -1 & 3 \\ \end{matrix} \right]$$
再消去 $y$：

② $\times3+$ ③得 $8z=24\quad ③$

即

$$\left[ \begin{matrix} x & y & z & w \\ 1 & 1 & 1 & 6 \\ 0 & -1 & 3 & 7 \\ 0 & 0 & 8 & 24 \\ \end{matrix} \right]$$

这时候我们就解出了，$z$然后回代到②解出$y$，再回代到①解出 x x <script type="math/tex" id="MathJax-Element-17">x</script>，这个方程组就解完了。

特殊情况

解方程的时候有时会遇到无解/多解的情况。

当方程组无解时，表现为消元时有一行所有系数都为0，但常数项不为0。

当方程组多解时，表现为消元时有一行所有系数都为0，且常数项也为0。有几行这种情况就说明有几个元的值可以随便取。

代码

以洛谷P3389为例：

因为在实际求解过程中会有实数存在，所以我们每次消元的时候选择系数最大的那一个，这样可以减少精度误差。

这道题没有无解的情况。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 105
#define eps 1e-8
using namespace std;
typedef double DB;
int n;
DB a[N][N],ans[N];
#define abs(x) ((x)<0?-(x):(x))
int main(){
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]);
    for (int i=1,x=1;i<=n;x=++i){//消元
        for (int j=i+1;j<=n;j++)
            if (abs(a[j][i])>abs(a[x][i])) x=j;
        if (abs(a[x][i])<eps) return puts("No Solution"),0;
        if (i!=x) swap(a[i],a[x]); DB d=a[i][i];
        for (int j=i;j<=n+1;j++) a[i][j]/=d;
        for (int j=i+1;j<=n;j++){
            d=a[j][i];
            for (int k=i;k<=n+1;k++)
                a[j][k]-=a[i][k]*d;
        }
    }
    for (int i=n;i;i--){//回代
        ans[i]=a[i][n+1];
        for (int j=i+1;j<=n;j++)
            ans[i]-=a[i][j]*ans[j];
    }
    for (int i=1;i<=n;i++)
        printf("%.2f\n",ans[i]);
    return 0;
}