高斯消元法模板

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

求解的过程中，在消去某个未知数时，打算保留该未知数的式子的对应未知数系数可能为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 105
double 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]
=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]
;
}
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]
=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]
;
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;
}