高斯消元法模板

来源：互联网发布：樱井知香喷泉gif 编辑：程序博客网时间：2024/06/06 18:36

高斯消元法可以用来求解线性方程组，在线性代数中学过，如果不了解，可参见高斯消元法-百度百科或高斯消元法-维基百科

求解的过程中，在消去某个未知数时，打算保留该未知数的式子的对应未知数系数可能为0，在这种情况下，只要调整方程的顺序，使得对应的系数不为0就可以了(另外，为了减小误差，应该选择要消去的未知数系数的绝对值尽可能大的方程，该方法被称为列主元高斯消元法)。

时间复杂度为O(n^3)。

理解了高斯消元法的原理，代码很容易看懂，代码里在关键地方有注释……

代码如下：

1.二维数组

#include <iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>using namespace std;const double EPS=1e-8;#define maxn 105double A[maxn][maxn];double B[maxn][maxn+1];double x[maxn];double b[maxn];int n;//求解Ax=b，其中A是方阵void gauss_jordan(){    for(int i=0;i<n;i++){        for(int j=0;j<n;j++){            B[i][j]=A[i][j];        }    }    for(int i=0;i<n;i++){        B[i][n]=b[i];    }    for(int i=0;i<n;i++){        //把正在处理的未知数系数的绝对值最大的式子换到第i行        int pivot=i;        for(int j=i;j<n;j++){            if(abs(B[j][i]>abs(B[pivot][i]))) pivot=j;        }        for(int j=0;j<=n;j++){            swap(B[i][j],B[pivot][j]);        }        //无解或者无穷多解        if(abs(B[i][i])<EPS) return;        //把正在处理的未知数的系数变成1        for(int j=i+1;j<=n;j++){            B[i][j]/=B[i][i];        }        for(int j=0;j<n;j++){            if(i!=j){                //从第j个式子中消去第i个未知数                for(int k=i+1;k<=n;k++){                    B[j][k]-=B[j][i]*B[i][k];                }            }        }    }    for(int i=0;i<n;i++){        x[i]=B[i][n];    }    return;}int main(){    cin>>n;    for(int i=0;i<n;i++){        for(int j=0;j<n;j++){            cin>>A[i][j];        }    }    for(int i=0;i<n;i++){        cin>>b[i];    }    gauss_jordan();    for(int i=0;i<n;i++){        cout<<x[i]<<" ";    }    printf("\n");    return 0;}

2.vector数组

#include <iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>#include<vector>using namespace std;const double EPS= 1e-8 ;typedef vector<double> vec;typedef vector<vec> mat;//求解Ax=b，其中A是方阵//当方程组无解或有无穷多解时，返回一个长度为0的数组vec gauss_jordan(const mat& A,const vec& b){    int n=A.size();    mat B(n,vec(n+1));    //把A复制给B    for (int i=0;i<n;i++)        for (int j=0;j<n;j++) B[i][j]=A[i][j];    //把正在处理的未知数系数的绝对值最大的式子换到第i行    for (int i=0;i<n;i++) B[i][n]=b[i];    for (int i=0;i<n;i++){        int pivot=i;        for (int j=i;j<n;j++){            if (abs(B[j][i])>abs(B[pivot][i])) pivot=j;        }        swap(B[i],B[pivot]);        //无解或有无穷多解        if (abs(B[i][i])<EPS) return vec();        //把正在处理的未知数系数变为1        for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i];        for (int j=0;j<n;j++){            if (i!=j){                //从第j个式子中消去第i个未知数                for (int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];            }        }    }    vec x(n);    //存放在右边的b就是答案    for (int i=0;i<n;i++) x[i]=B[i][n];    return x;}int main(){    int n;    cin>>n;    mat A(n,vec(n));    vec b(n);    for(int i=0;i<n;i++){        for(int j=0;j<n;j++){            cin>>A[i][j];        }    }    for(int i=0;i<n;i++){        cin>>b[i];    }    vec x=gauss_jordan(A,b);    for(int i=0;i<n;i++){        cout<<x[i]<<" ";    }    cout<<endl;    return 0;}

0 0