淺談高斯消元

Jack_Pei發表於2019-05-10

高斯消元。。。當初以為自己學會了。。。後來。。。

話說這個東西好像最早出現於《九章算術》誒（古代人就是強）

廢話不說，進入正題。。。

前置知識

高斯消元法是解線性方程組的方法之一 

首先，線性方程組瞭解一下：

　　可認為線性方程組就是一次方程組。如圖：

 

如果存在常數c1,c2,c3,...,cn代替x1,x2,x3,...,xn,使上圖的每個方程成立，則稱（c1,c2,c3,...,cn）為方程組的一個解. 解方程就是求其全部解

顯然對線性方程組進行以下三種變換,所得方程組解集不變(這三種變換被稱為同解變換)

　　(1) 對換兩個方程（換法變換）（其實就是交換兩個方程的位置，上下交換) ;

　　(2) 用非零數 c 乘以某一個方程（倍法變換）;

　　(3) 把某一個方程的 k 倍加到另一個方程上去( 消法變換) .

的確顯然吧...

 

再來了解一下矩陣(=_=)：

由 s*n 個數 aij ( i = 1, 2, … , s; j = 1, 2, … , n) 排成的矩形表(可理解為長是n，寬是s的二維陣列。。)，稱為 s 行 n 列矩陣 .aij稱為矩陣的( i , j) 元 .通常用大寫字母或 ( aij )sn 表示矩陣 .如果 s = n（長等於寬）稱A是 n 階方陣或 n 階矩陣 .

全零矩陣記做Osn,或O；

n階單位矩陣記做Inn,或I（單位矩陣是個方陣，從左上角到右下角的對角線（稱為主對角線）上的元素均為1。除此以外全都為0。）

 

對於矩陣(aij)sn、(bij)ns，若aij=bji(1<=i<=s,1<=j<=n)則稱B為A的轉置矩陣，記做B=AT

 

另：係數矩陣和增廣矩陣

把線性方程組的係數aij列入矩陣，稱為係數矩陣。

如圖：

而在係數矩陣右側加入一列常數項，稱其為增廣矩陣。

如圖：

（好吧增廣矩陣就是類似把未知數扣掉了再扔到矩陣裡。。）

對矩陣的行(列)做 線性方程組 的三種同解變換（上面的），稱為矩陣的初等變換。

 

高斯消元法

對方程組的增廣矩陣做初等變換，將係數部分化為對角線，上三角形，或階梯形

思路：找到一個不為零的係數，設其所在行為pos，列為i，然後拿這個數去消掉其它行，但同列的數；但是消元的時候要有一定順序

步驟：

先把一個個方程扔到矩陣中，構建增廣矩陣

然後從第i行（i從1到n，按順序）開始，在i+1行到n中的第i列（上面的方程已經符合上三角性質了）選取一個不為零的係數，記其所在行為pos，

交換第pos行與第i行，然後pos行變成了第i行（注：後面的第i行為交換完後的，即交換前的pos行）

至於為什麼要找最大的，這樣有利於提高精度。。（學長說的。。）

於是拿第i行的第i項去消i+1行到n行的第i項

這樣可以消成上三角矩陣；

然後從最後一行開始，向上消去不需要的係數（最後一行只有一個係數，可以直接計算答案）

程式碼（Luogu 模板）：

#include<cstdio>
#include<iostream>
#define R register int
using namespace std;
const double eps=1E-11;
inline int g() {
    R ret=0,fix=1; register char ch; while(!isdigit(ch=getchar())) fix=ch=='-'?-1:fix;
    do ret=ret*10+(ch^48); while(isdigit(ch=getchar())); return ret*fix;
}
int n;
double a[110][110];
inline bool ck0(double x) {return x<eps&&x>-eps;}
inline void Gauss() {
    for(R i=1,pos;i<=n;++i) {
        pos=i; for(R j=i+1;j<=n;++j) if(ck0(a[pos][i])||a[pos][i]<a[j][i]) pos=j;
        if(pos!=i) for(R j=1;j<=n+1;++j) swap(a[i][j],a[pos][j]);
        if(ck0(a[i][i])) {printf("No Solution\n"); return ;}
        for(R j=i+1;j<=n;++j) {
            register double tmp=a[j][i]/a[i][i];
            for(R k=1;k<=n+1;++k) a[j][k]-=tmp*a[i][k];
        }
    } for(R i=n;i>=1;--i) {
        for(R j=n;j>i;--j) a[i][n+1]-=a[i][j]*a[j][n+1];//把這行係數代表的答案消掉
        a[i][n+1]/=a[i][i];
    } for(R i=1;i<=n;++i) printf("%.2lf\n",a[i][n+1]);
}
signed main() {
    n=g(); for(R i=1;i<=n;++i) for(R j=1;j<=n+1;++j) a[i][j]=g();
    Gauss();
}

 

(例題)

Luogu P3389 【模板】高斯消元法

Luogu P4035 [JSOI2008]球形空間產生器

POJ1830   

etc. 

2019.05.10

 

 

 

 
 

相關文章