高斯消元

$转载自 https://www.luogu.com.cn/blog/pks-LOVING/p3389-mu-ban-gao-si-xiao-yuan-fa$

$e m m m m 这个消元方式其实严格来说是可行性算法\dots\dots而不是优化性算法\dots\dots不过话说由于我太蒟蒻了，所以并不知道什么更优的算法（#滑稽）$

嗯，其实这个算法是 $O(n^3)O(n3) 的算法，需要一些矩阵及行列式的知识……那么由本蒟蒻来记录一下这个算法吧！$

那么假设有一个线性方程组是长这样的：

$\5x & + & y &+& 6z &=&25 \2x & + & 3y &+& 4z &=&20end{Bmatrix}⎩⎨⎧3x5x2x+++2yy3y+++z6z4z===102520⎭⎬⎫$

初中老师说过：我们可以加减消元或者代入消元,但是我们需要在程序里实现的时候，需要一种有规律可循的算法。所以我们选择加减消元，但用代入消元回带。

整体思路就是我们可以先在某一个式子里，用这个式子的

比如刚才这个方程，我们用第二个式子去消1、3式里的

$\5x & + & y &+& 6z &=&25 \0 imes x & + & frac{13}{5}y &+& frac{8}{5}z &=&10end{Bmatrix}⎩⎨⎧0×x5x0×x+++57yy513y+++(−513z)6z58z===−52510⎭⎬⎫$

整理之后再用第三个式子里的

$\ frac{13}{5}y &+& frac{8}{5}z &=&10end{Bmatrix}{0×y513y++(−65225z)58z==−1313510}$

那么我们发现在

带回三式里解出

再将

好像这个方法再数学逻辑上讲是特别麻烦的，但是却是一个通用性强的方法

那么放在程序实现上来讲，我们可以先用一个

$\5 & 1 & 6 &|& 25 \2 & 3 & 4 &|&20end{Bmatrix}⎩⎨⎧352213164∣∣∣102520⎭⎬⎫$

左边记录系数，右边记录每个式子的结果。

那么首先我们需要用第一列中（所有的

由于最多也只能用

假设我们现在在处理第 $k_1 k_2 k_3 k_ 4k1k2k3k4 …… kmkm ,那么考虑，在选完 k_iki 之后，下面我们要进行的是把 k_iki 消成 11 。那么此时对于第 ii 行的其他的系数以及结果我们都要除以 k_iki 。$

之后呢？之后我们要进行的操作是用这个式子来消掉其他式子里的该未知数啊 $k_iki 的式子,则对于每个式子 jj 而言，每个系数减去当前系数的倍数，这个倍数应该为 k_jkj ）$

那么这样看来，对于当次用来消元的式子的每个系数 $q_{i1}q_{i2}q_{i3}q_{i4}qi1qi2qi3qi4 …… q_{iw}qiw (假设当前元的系数是 q_{i1}qi1 )而言，对于每一个其他式子的该项系数 q_{jw}qjw ,都需要让 q_{jw}qjw 变成$

$q_{jw}-frac{q_{j1}}{q_{i1}} imes q_{iw}qjw−qi1qj1×qiw$

那么我们观察这个式子， $q_{i1}qi1 越大， frac{q_{j1}}{q_{i1}}qi1qj1 期望越小，那么我们考虑，这个值越小，我们就约可以把它看作一个“基本单位”。从而我们就使得减出来的值失精程度越低，最后即可保证数据是从期望上来讲最精确。$

嗯，讲的很麻烦，大家挑重点看吧（或者只看最后一个自然段）

(逃

在置为

代码长这个样子：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
double map[111][111];
double ans[111];
double eps=1e-7;
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&map[i][j]);
    for(int i=1;i<=n;i++){
        int r=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(map[r][i])<fabs(map[j][i]))
                r=j;//find_the_biggest_number_of_the_first_column（at present) 
        if(fabs(map[r][i])<eps){
            printf("No Solution");
            return 0;
        }
        if(i!=r)swap(map[i],map[r]);//对换一行或一列,属于找最大当前系数的其中一步。（这样就可以只处理当前行的系数啦！） 
        double div=map[i][i];
        for(int j=i;j<=n+1;j++)
            map[i][j]/=div;
        for(int j=i+1;j<=n;j++){
            div=map[j][i];
            for(int k=i;k<=n+1;k++)
                map[j][k]-=map[i][k]*div;
        }
    }
    ans[n]=map[n][n+1];
    for(int i=n-1;i>=1;i--){
        ans[i]=map[i][n+1];
        for(int j=i+1;j<=n;j++)
            ans[i]-=(map[i][j]*ans[j]);
    }//回带操作
    for(int i=1;i<=n;i++)
        printf("%.2lf
",ans[i]);
}