天天看点

高斯消元(洛谷P3389)简介实现代码

简介

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

实现

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

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

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

举个例子

比如下面这个方程组:

⎧⎩⎨x+y+z2x+3y−zx−2y+2z=6①=5②=3③(1) (1) { x + y + z = 6 ① 2 x + 3 y − z = 5 ② x − 2 y + 2 z = 3 ③

我们用矩阵来表示:

⎡⎣⎢⎢⎢x121y13−2z1−11w650⎤⎦⎥⎥⎥ [ x y z w 1 1 1 6 2 3 − 1 5 1 − 2 1 0 ]

现在我们消去 x x :

① ×2−×2− ②得  −y+3z=7②   − y + 3 z = 7 ②

① − − ③得 3y−z=3③3y−z=3③

⎡⎣⎢⎢⎢x100y1−13z13−1w673⎤⎦⎥⎥⎥ [ x y z w 1 1 1 6 0 − 1 3 7 0 3 − 1 3 ]

再消去 y y :

② ×3+×3+ ③得 8z=24③ 8 z = 24 ③

⎡⎣⎢⎢⎢x100y1−10z138w6724⎤⎦⎥⎥⎥ [ x y z w 1 1 1 6 0 − 1 3 7 0 0 8 24 ]

这时候我们就解出了, z z 然后回代到②解出yy,再回代到①解出 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=;i<=n;i++)
        for (int j=;j<=n+;j++)
            scanf("%lf",&a[i][j]);
    for (int i=,x=;i<=n;x=++i){//消元
        for (int j=i+;j<=n;j++)
            if (abs(a[j][i])>abs(a[x][i])) x=j;
        if (abs(a[x][i])<eps) return puts("No Solution"),;
        if (i!=x) swap(a[i],a[x]); DB d=a[i][i];
        for (int j=i;j<=n+;j++) a[i][j]/=d;
        for (int j=i+;j<=n;j++){
            d=a[j][i];
            for (int k=i;k<=n+;k++)
                a[j][k]-=a[i][k]*d;
        }
    }
    for (int i=n;i;i--){//回代
        ans[i]=a[i][n+];
        for (int j=i+;j<=n;j++)
            ans[i]-=a[i][j]*ans[j];
    }
    for (int i=;i<=n;i++)
        printf("%.2f\n",ans[i]);
    return ;
}