淺談高斯消元的實現和簡單應用

1、高斯消元的原理

對於n元的m個線性方程組成的方程組，咱們將其以矩陣的形式記錄下來：

a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn

而後進行初等行列變換，嘗試構造出一個上三角矩陣，逐步使係數不爲零的項減小；

等最後只剩下一個係數不爲零時，進行回代，逐步求出已知解。（詳解過程諮詢小學老師）

2、高斯消元的實現

老老實實的回代代碼參見其餘人的博客，這裏介紹一種比較毒瘤的不回代暴力消元法：

1.Process

對於每一個方程，按照必定的規則（後話）挑選一個主元，記錄該主元對應第幾個方程，而後用初等行列變換消去其餘全部該元的係數；

最後咱們獲得的是一個每行只有一個數不爲零，每列只有一個數不爲零的鬼畜矩陣（本身腦補）

此時令ans向量對應的數字出去該行的非零係數，即爲對應該元的解。

2.Code

設a數組爲已經記錄係數的數組（格式見上方），那麼a應該是n行n+1列的，最後一列表明方程的常數項；

設w數組記錄每一個方程的主元是第幾項，v數組記錄答案；

當多解時輸出「Multiple solutions」，無解時輸出」No solution」；

double a[max_n][max_n+1],v[max_n];

int w[max_n]; 

void gauss(){
    double eps=1e-6;
    for(int i=1;i<=n;++i){    //Enumerate the equation;
        int p=0;                //Record the position of the largest number; 
        double mx=0;        //Recording the largest number;
        for(int j=1;j<=n;++j)
            if(fabs(a[i][j])-eps>mx){
                mx=fabs(a[i][j]);p=j;    //fabs() returns the absolute value of float; 
            }
        if(!p){
            if(fabs(a[i][n+1]<eps))printf("Multiple solutions");
            else printf("No solution");
            return; 
        }
        w[i]=p;
        for(int j=1;j<=n;++j)
            if(i!=j){       //other equations
                double t=a[j][p]/a[i][p];
                for(int k=1;k<=n+1;++k)    //n+1 is important
                    a[j][k]-=a[i][k]*t;
            }
    }
    for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
}

3.notice

（1）精度的設置

衆所周知浮點數是有精度丟失的，在高斯消元中，精度的選擇要依題目而定，精度太低會致使係數較小的數被誤認爲係數爲零，而精度太高也有可能會致使偏差大於精度，致使本應該係數消爲0的項誤認爲係數不爲零，因此精度的選擇是很哲學的問題，要注意。

推薦範圍：1e-4到1e-10

（2）主元的選取原則

接着（1）說，咱們知道，用浮點數是有精度丟失的，那麼讓一個較大的數除以一個極小的數偏差天然大的可想而知，因此咱們想獲得在精度容許的條件下係數最大的主元，因此對於每一個方程，咱們都應該選擇最大系數的元做爲主元。

（3）在模2意義下的高斯消元

使用bitset優化運行時間，詳見相關應用中第三個例題的講解；

3、相關應用

這裏給出高斯消元的幾道基礎題目，難度適合初學者。

1.[Luogu P3389]【模板】高斯消元

Description

給定一個線性方程組，對其求解

輸入格式：
第一行，一個正整數 n
第二至 n+1行，每行 n+1個整數，爲 a1,a2⋯an和 b，表明一組方程。

輸出格式：
共n行，每行一個數，第 i行爲 xi(保留2位小數）
若是不存在惟一解，在第一行輸出"No Solution".

Solution

如上所述。

Code

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;

const int max_n=110;
double a[max_n][max_n+1],v[max_n];
int n,w[max_n]; 

inline int rd(){
    int x=0;
    bool f=0;
    char c=getchar();
    while(!isdigit(c)){
        if(c=='-') f=1;
        c=getchar();
    }
    while(isdigit(c)){
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return f?-x:x;
}

void gauss(){
    double eps=1e-6;
    for(int i=1;i<=n;++i){//enumerate the equation;
        int p=0;          //Record the position of the largest number; 
        double mx=0;      //Recording the largest number;
        for(int j=1;j<=n;++j)
            if(fabs(a[i][j])-eps>mx){
                mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; 
            }
        if(!p){
            printf("No Solution");
            return; 
        }
        w[i]=p;
        for(int j=1;j<=n;++j)
            if(i!=j){       //other equations
                double t=a[j][p]/a[i][p];
                for(int k=1;k<=n+1;++k)//n+1 is important
                    a[j][k]-=a[i][k]*t;
            }
    }
    for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
    for(int i=1;i<=n;++i) printf("%.2lf\n",v[i]);
}

int main(){
    n=rd();
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n+1;++j)
            a[i][j]=rd();
    gauss();
    return 0; 
}

2.[BZOJ 1013][JSOI 2008] 球形空間產生器sphere

詳解參考個人隨筆：http://www.cnblogs.com/COLIN-LIGHTNING/p/8982341.html