高斯消元法

根据数值分析的高斯消元算法，可写出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