高斯消元法

根据数值分析的高斯消元算法，可写出C++的实现代码，如下所示：

#include<bits/stdc++.h>using namespace std;const int maxn=105;double a[maxn][maxn],m[maxn][maxn],b[maxn],x[maxn],det=1.0;int n;void input_data(){freopen("data.txt","r",stdin);cin>>n;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}}for(int i=1;i<=n;i++) {cin>>b[i];}}void print(){for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){printf("%.4lf ",a[i][j]);}printf("%.4lf\n",b[i]);}}void gauss(){for(int k=1;k<n;k++)//枚举每一行 {printf("===========\n");//我是分割线 print();//打印出消元过程，系数矩阵的变化 det*=a[k][k];for(int i=k+1;i<=n;i++)//利用第k行消第i行 {m[i][k]=a[i][k]/a[k][k];a[i][k]=0;for(int j=k+1;j<=n;j++)// 消第i行的第k列 {a[i][j]=a[i][j]-m[i][k]*a[k][j];}b[i]=b[i]-m[i][k]*b[k];}}x[n]=b[n]/a[n][n];//得到解x[n] for(int i=n-1;i>=1;i--)//依次求出解 {double sum=0.0;for(int j=i+1;j<=n;j++){sum+=a[i][j]*x[j];}x[i]=(b[i]-sum)/a[i][i];}printf("===========\n");print();det*=a[n][n];}void output_ans(){cout<<"ans is :"<<endl;for(int i=1;i<=n;i++) cout<<x[i]<<" ";cout<<endl;}void output_LU()//得到LU矩阵 {cout<<"Lower triangular matrix is :"<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){if(i==j) cout<<1<<" ";else if(i>j) cout<<m[i][j]<<" ";else cout<<0<<" ";}cout<<endl;}cout<<"Upper triangular matrix"<<endl;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){if(i<=j) cout<<a[i][j]<<" ";else cout<<0<<" ";}cout<<endl;}}void output_det()//得到行列式的值 {cout<<"the det value of matrix is :"<<endl;cout<<det<<endl;}int main(){//from text input data input_data();gauss();output_ans();return 0;} 

进一步，优化高斯消元法，可以得到列主元消元法

#include<bits/stdc++.h>using namespace std;const int maxn=105;double a[maxn][maxn],m[maxn][maxn],b[maxn],x[maxn],det=1.0;int n;void input_data(){freopen("data2.txt","r",stdin);cin>>n;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}}for(int i=1;i<=n;i++) cin>>b[i];}void print(){for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){printf("%.4lf ",a[i][j]);}printf("%.4lf\n",b[i]);}}void column_main_element(){for(int k=1;k<n;k++){printf("===========\n");print();int pos=k;for(int i=k+1;i<=n;i++)if(a[i][k]>a[pos][k]) pos=i;for(int i=k;i<=n;i++)swap(a[k][i],a[pos][i]);swap(b[k],b[pos]);det*=a[k][k];if(pos!=k) det*=-1;if(fabs(a[k][k])<1e-6){det=0.0;return ;}printf("==========\n");print();for(int i=k+1;i<=n;i++){m[i][k]=a[i][k]/a[k][k];a[i][k]=0;for(int j=k+1;j<=n;j++){a[i][j]=a[i][j]-m[i][k]*a[k][j];}b[i]=b[i]-m[i][k]*b[k];}}x[n]=b[n]/a[n][n];for(int i=n-1;i>=1;i--){double sum=0.0;for(int j=i+1;j<=n;j++){sum+=a[i][j]*x[j];}x[i]=(b[i]-sum)/a[i][i];}det*=a[n][n];printf("===========\n");print();}void output_ans(){for(int i=1;i<=n;i++) cout<<x[i]<<" ";cout<<endl;}int main(){input_data();column_main_element();output_ans();return 0;}

最后，利用高斯消元法，可以求逆矩阵。即gauss-jordan算法

#include<bits/stdc++.h>using namespace std;const int maxn=105;double a[maxn][maxn],m[maxn][maxn],b[maxn],x[maxn],det=1.0;int n;void input_data(){freopen("data4.txt","r",stdin);cin>>n;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}a[i][n+i]=1;}for(int i=1;i<=n;i++) cin>>b[i];}void print(){for(int i=1;i<=n;i++){for(int j=1;j<=2*n;j++){printf("%.4lf ",a[i][j]);}printf("%.4lf\n",b[i]);}}void column_main_element(){for(int k=1;k<n;k++){//printf("===========\n");//print();int pos=k;for(int i=k+1;i<=n;i++)if(a[i][k]>a[pos][k]) pos=i;for(int i=k;i<=2*n;i++)swap(a[k][i],a[pos][i]);swap(b[k],b[pos]);det*=a[k][k];if(pos!=k) det*=-1;if(fabs(a[k][k])<1e-6){det=0.0;return ;}m[k][k]=1.0/a[k][k];for(int i=k+1;i<=n;i++){m[i][k]=a[i][k]/a[k][k];a[i][k]=0;for(int j=k+1;j<=2*n;j++){a[i][j]=a[i][j]-m[i][k]*a[k][j];}b[i]=b[i]-m[i][k]*b[k];}for(int j=k;j<=2*n;j++){a[k][j]=a[k][j]*m[k][k];}b[k]=b[k]*m[k][k];printf("===========\n");print();}m[n][n]=1.0/a[n][n];a[n][n]=a[n][n]*m[n][n];b[n]=b[n]*m[n][n];x[n]=b[n]/a[n][n];for(int p=n+1;p<=2*n;p++) {a[n][p]*=m[n][n];}//for(int i=n-1;i>=1;i--)//{//double sum=0.0;//for(int j=i+1;j<=n;j++)//{//sum+=a[i][j]*x[j];//}//x[i]=(b[i]-sum)/a[i][i];//}for(int i=n-1;i>=1;i--){for(int j=i;j>=1;j--){double tmp=-a[j][i+1];b[j]+=tmp*b[i+1];a[j][i+1]=0;for(int p=n+1;p<=2*n;p++){a[j][p]+=tmp*a[i+1][p]; } } x[i]=b[i];}det*=a[n][n];printf("===========\n");print();}void output_ans(){for(int i=1;i<=n;i++) cout<<x[i]<<" ";cout<<endl;}int main(){input_data();printf("==========\n");print();column_main_element();output_ans();return 0;}

实验数据1，解线性方程组，C++代码的输入：

3
1 1 1
0 4 -1
2 -2 1
6 5 1

其含义是

 x11   +   x12  +  x13=6

            4*x22  -  x23=5

2*x31 - 2*x32 +  x33=1

其得到的输出为：

===========
1.0000 1.0000 1.0000 6.0000
0.0000 4.0000 -1.0000 5.0000
2.0000 -2.0000 1.0000 1.0000
===========
1.0000 1.0000 1.0000 6.0000
0.0000 4.0000 -1.0000 5.0000
0.0000 -4.0000 -1.0000 -11.0000
===========
1.0000 1.0000 1.0000 6.0000
0.0000 4.0000 -1.0000 5.0000
0.0000 0.0000 -2.0000 -6.0000
ans is :
1 2 3