高斯消元和高斯约旦消元 Gauss（-Jordan） Elimination

高斯消元法，是线性代数中的一个算法，可用来求解线性方程组，并可以求出矩阵的秩，以及求出可逆方阵的逆矩阵。

在讲算法前先介绍些概念

矩阵的初等变换

矩阵的初等变换又分为矩阵的初等行变换和矩阵的初等列变换。矩阵的初等行变换和初等列变换统称为初等变换。另外：分块矩阵也可以定义初等变换。

等价

定义：如果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;
}

 例题

[SDOI2007] 线性方程组

代码

#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
",i+1);
        else printf("x%d=%.2lf
",i+1,ans[i]);
    }
    return 0;
}