高斯消元法,是線性代數中的一個算法,可用來求解線性方程組,並可以求出矩陣的秩,以及求出可逆方陣的逆矩陣。
在講算法前先介紹些概念
定義:如果B可以由A經過一系列初等變換得到,則稱矩陣A與B稱為等價
初等行變換
定義:所謂數域P上矩陣的初等行變換是指下列3種變換:
1)以P中一個非零的數乘矩陣的某一行
2)把矩陣的某一行的c倍加到另一行,這里c是P中的任意一個數
3)互換矩陣中兩行的位置
一般來說,一個矩陣經過初等行變換后就變成了另一個矩陣,當矩陣A經過初等行變換變成矩陣B時,一般寫作
可以證明:任意一個矩陣經過一系列初等行變換總能變成
階梯型矩陣。(方陣的話就是上三角矩陣)
初等列變換
同樣地,定義初等列變換,即:
1)以P中一個非零的數乘矩陣的某一列
2)把矩陣的某一列的c倍加到另一列,這里c是P中的任意一個數
3)互換矩陣中兩列的位置
初等矩陣
定義:由單位矩陣E經過一系列初等變換得到的矩陣稱為初等矩陣。
引理:對一個
矩陣A作一初等行變換,就相當於在A的左邊乘上相應的
的初等矩陣;
對A作一初等列變換,就相當於在A的右邊乘上相應的
的初等矩陣。
高斯消元的思想——消元法(加減消元 & 代入消元)
消元法是將方程組中的一方程的未知數用含有另一未知數的代數式表示,並將其代人到另一方程中,這就消去了一未知數,得到一解;或將方程組中的一方程倍乘某個常數加到另外一方程中去,也可達到消去一未知數的目的。消元法主要用於二元一次方程組的求解。
步驟
首先要把方程轉化為增廣矩陣矩陣,比如
將上面的方程轉為下面的增廣矩陣
其中第4列是常數項,其它三列是方程的系數。
然后開始進行消元
高斯消元
首先找到一個第一列的數不為0的行(一般找第一列的數最大的行)(如果都為0就跳過當前步)
然后用它的第一列的數將下面行當前列的值化為0,變換過的初等矩陣與原矩陣等價,化為方程后依然成立
本矩陣第一列的數最大的為第三行就把第三行與第一行交換

然后下面行的當前列消去,如圖

除了最后一列外,每一列都如此,最后得到上三角矩陣如圖

這樣我們很容易算出x3的值,再用x3的值算出x2,x1的值
int gauss(int n){//a為增廣矩陣 int r,w; for(int i=0;w<n&&i<n;i++,w++){//進行到第i列,第w行 int r=w; for(int j=w+1;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到當前列絕對值最大的行 if(fabs(a[r][i])<eps){w--;continue;}//當前列值都為0,跨過當前步 if(w!=r)for(int j=i;j<=n;j++)swap(a[r][j],a[w][j]);//交換當前列絕對值最大的行和沒計算過的第一行 for(int k=w+1;k<n;k++){//消去當前列(只消下面行的) double pro=a[k][i]/a[w][i]; for(int j=i;j<=n;j++) a[k][j]-=pro*a[w][j]; } } return w; }
高斯約旦消元
同高斯消元大體差不多,但消元時上面計算過的行也要消去當前列,最后得到的是對角矩陣而不是上三角矩陣。
第一次和高斯消元相同

第二次要順帶把第一行第二列消去

最后將最后一行上面所有行的倒數第二列消去

int gauss_jordan(int n){//a為增廣矩陣 int r,w=0; for(int i=0;i<n&&w<n;w++,i++){//進行到第i列,第w行 int r=w; for(int j=w+1;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到當前列絕對值最大的行 if(fabs(a[r][i])<eps){w--;continue;}//當前列值都為0,跨過當前步 if(r!=w)for(int j=0;j<=n;j++)swap(a[r][j],a[w][j]);//交換當前列絕對值最大的行和沒計算過的第一行 for(int k=0;k<n;k++){//消去當前列(除本行外) if(k!=w) for(int j=n;j>=w;j--)a[k][j]-=a[k][i]/a[w][i]*a[w][j]; } } return w; }
例題
代碼
#include<cstdio> #include<cmath> #include<algorithm> #define eps 1e-8 using namespace std; double a[55][55],ans[55]; int d; int gauss_jordan(int n){//a為增廣矩陣 int r,w=0; for(int i=0;i<n&&w<n;w++,i++){//進行到第i列,第w行 int r=w; for(int j=w+1;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到當前列絕對值最大的行 if(fabs(a[r][i])<eps){w--;continue;}//當前列值都為0,跨過當前步 if(r!=w)for(int j=0;j<=n;j++)swap(a[r][j],a[w][j]);//交換當前列絕對值最大的行和沒計算過的第一行 for(int k=0;k<n;k++){//消去當前列(除本行外) if(k!=w) for(int j=n;j>=w;j--)a[k][j]-=a[k][i]/a[w][i]*a[w][j]; } } return w; } int gauss(int n){//a為增廣矩陣 int r,w; for(int i=0;w<n&&i<n;i++,w++){//進行到第i列,第w行 int r=w; for(int j=w+1;j<n;j++)if(fabs(a[j][i])>fabs(a[r][i]))r=j;//找到當前列絕對值最大的行 if(fabs(a[r][i])<eps){w--;continue;}//當前列值都為0,跨過當前步 if(w!=r)for(int j=i;j<=n;j++)swap(a[r][j],a[w][j]);//交換當前列絕對值最大的行和沒計算過的第一行 for(int k=w+1;k<n;k++){//消去當前列(只消下面行的) double pro=a[k][i]/a[w][i]; for(int j=i;j<=n;j++) a[k][j]-=pro*a[w][j]; } } return w; } int main(){ freopen("gaess.in","r",stdin); // freopen("gaess.out","w",stdout); int n;scanf("%d",&n); for(int i=0;i<n;i++){ for(int j=0;j<=n;j++){ scanf("%lf",a[i]+j); } } #if 0 //guass d=gauss(n);d--; for(int j=d;j<n;j++){ if(fabs(a[j][n-1])<eps&&fabs(a[j][n])>eps){printf("-1");return 0;}//最后一個方程未知數最少,方程=右邊不為0=左邊為0,則無解 } if(d<n-1){printf("0");return 0;}//一個方程解一個未知數,有效方程少於n個,則多個解 for(int i=d;i>=0;i--){ for(int k=i+1;k<n;k++)a[i][n]-=a[i][k]*ans[k]; ans[i]=a[i][n]/a[i][i]; } #endif #if 1 //guass_jordan d=gauss_jordan(n);d--; for(int i=0;i<n;i++){//有一個方程=右邊不為0=左邊為0,則無解 bool d=1; for(int j=i;j<n;j++)d&=(fabs(a[i][j])<eps); if(d&&fabs(a[i][n])>eps){ printf("-1");return 0; } } for(int i=0;i<n;i++){//消元后有變量在多個方程中出現,則有多個解 int max1=0; for(int j=i;j<n;j++)if(fabs(a[i][j])>eps)max1++; if(max1>1){printf("0");return 0;} } for(int i=0;i<n;i++){ ans[i]=a[i][n]/a[i][i]; } #endif for(int i=0;i<n;i++){ if(fabs(ans[i])<eps)printf("x%d=0\n",i+1); else printf("x%d=%.2lf\n",i+1,ans[i]); } return 0; }
