高斯列主消元詳解及模板

採用高斯先列主元消元法求解線性方程組AX=b

方法說明（以4階為例）：

（1）第1步消元——在增廣矩陣（A，b）第一列中找到絕對值最大的元素，將其所在行與第一行交換，再對（A，b）

做初等行變換使原方程組轉化為如下形式：注：“*”代表非0。

*x1+*x2+*x3+*x4=y1;
 0 +*x2+*x3+*x4=y2;
 0 +*x2+*x3+*x4=y3;
 0 +*x2+*x3+*x4=y4;

（2）第2步消元——在增廣矩陣（A，b）中的第二列中（從第二行開始）找到絕對值最大的元素，將其所在行與第二行交換，再對（A，b）做初等行變換使原方程組轉化為：

*x1+*x2+*x3+*x4=y1;
 0 +*x2+*x3+*x4=y2;
 0 + 0 +*x3+*x4=y3;
 0 + 0 +*x3+*x4=y4;

（3）第3步消元——在增廣矩陣（A，b）中的第三列中（從第三行開始）找到絕對值最大的元素，將其所在行與第二行交換，再對（A，b）做初等行變換使原方程組轉化為：

*x1+*x2+*x3+*x4=y1;
 0 +*x2+*x3+*x4=y2;
 0 + 0 +*x3+*x4=y3;
 0 + 0 + 0 +*x4=y4;

（4）按x4 ; x3; ; x1 的順序回代求解出方程組的解

再回代求解即可。

下面是高斯消元模板：

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#include<cmath>
#define eps 0.0000001
using namespace std;
double x[102];
double g[102][102];
int flag;

void gauss(int n,int m)
{
    int row,col,i,j,k;
    for(row=1,col=1;row<n,col<m;row++,col++)
    {
        k=row;
        for(i=row+1;i<=n;i++)  //列主元
            if(fabs(g[i][col])>fabs(g[k][col]))
                k=i;
        if(k!=row)   //行交換
        {
            for(i=col; i<=m; i++)
                swap(g[k][i],g[row][i]);
        }
        if(fabs(g[row][col])<eps) //主元是0
        {
            flag=1;
            return;
        }
        for(i=row+1; i<=n; i++)  //主元不是0把下面的行第一個值全部變為0
        {
            if(fabs(g[i][col])<eps)
                continue;
            double t=g[i][col]/g[row][col];
            g[i][col]=0.0;
            for(j=col+1;j<=m;j++)
                g[i][j]-=t*g[row][j];
        }
    }
    if(fabs(g[n][n])<eps)
    {
        flag=1;
        return;
    }
    for(i=n;i>=1;i--)  //回代求解
    {
        x[i]=g[i][m];
        for(j=i+1;j<=n;j++)
            x[i]-=x[j]*g[i][j];
        x[i]/=g[i][i];
    }
}

int main()
{
    int n,m,i,j;
    while(~scanf("%d%d",&n,&m))
    {
        flag=0;
        for(i=1;i<=n;i++)
            for(j=1;j<=m;j++)
              scanf("%lf",&g[i][j]);
              
        gauss(n,m);
        if(flag)
        {
            cout<<"方程組無解或有無陣列解"<<endl;
            continue;
        }
        cout<<"方程組的解為:"<<endl;
        for(i=1;i<=n;i++)
            cout<<"x"<<i<<":   "<<x[i]<<endl;
    }
    return 0;
}

/*
2 3
2 3 5
1 4 6
2 3
2 3 5
0 0 6
2 3
2 3 5
0 4 6
*/

再推薦一首歌：

這一路走來來自楊宗緯