跳到主要内容

消元

{a1,1x1+a1,2x2+…+a1,nxn=a1,n+1a2,1x1+a2,2x2+…+a2,nxn=a2,n+1…an,1x1+an,2x2+…+an,nxn=an,n+1\begin{cases} a_{1,1} x_1 + a_{1,2} x_2 + \ldots + a_{1,n} x_n = a_{1,n+1} \\ a_{2,1} x_1 + a_{2,2} x_2 + \ldots + a_{2,n} x_n = a_{2,n+1} \\ \ldots \\ a_{n,1} x_1 + a_{n,2} x_2 + \ldots + a_{n,n} x_n = a_{n,n+1} \\ \end{cases}

参考资料​

算法对比​

算法名称时间复杂度代码长度需要回代区分无解和无数解
高斯 - 约旦消元法O(n3)O(n^3)约 2525 行否困难
高斯消元法O(n3)O(n^3)约 4040 行是容易

高斯 - 约旦消元法​

{a1,1a1,2…a1,n∣a1,n+1a2,1a2,2…a2,n∣a2,n+1…an,1an,2…an,n∣an,n+1→{10…0∣a1,n+1′01…0∣a2,n+1′…00…1∣an,n+1′\left\{\begin{matrix} a_{1,1} & a_{1,2} & \ldots & a_{1,n} \mid a_{1,n+1} \\ a_{2,1} & a_{2,2} & \ldots & a_{2,n} \mid a_{2,n+1} \\ \ldots & & & \\ a_{n,1} & a_{n,2} & \ldots & a_{n,n} \mid a_{n,n+1} \end{matrix}\right. \rightarrow \left\{\begin{matrix} 1 & 0 & \ldots & 0 \mid a'_{1,n+1} \\ 0 & 1 & \ldots & 0 \mid a'_{2,n+1} \\ \ldots & & & \\ 0 & 0 & \ldots & 1 \mid a'_{n,n+1} \end{matrix}\right. 
double a[N][N];
bool gauss(int n)
{
	for(int i=1;i<=n;i++)
	{
		int t=i;
		for(int j=i+1;j<=n;j++)
		{
			if(fabs(a[t][i])<fabs(a[j][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		double tmp=a[i][i];
		for(int j=1;j<=n+1;j++)a[i][j]/=tmp;
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			tmp=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;k++)
			{
				a[j][k]-=a[i][k]*tmp;
			}
		}
	}
	return 1;
}

高斯消元法​

{a1,1a1,2…a1,n∣a1,n+1a2,1a2,2…a2,n∣a2,n+1…an,1an,2…an,n∣an,n+1→{a1,1′a1,2′…a1,n∣a1,n+1′0a2,2…a2,n′∣a2,n+1′…00…an,n′∣an,n+1′\left\{\begin{matrix} a_{1,1} & a_{1,2} & \ldots & a_{1,n} \mid a_{1,n+1} \\ a_{2,1} & a_{2,2} & \ldots & a_{2,n} \mid a_{2,n+1} \\ \ldots \\ a_{n,1} & a_{n,2} & \ldots & a_{n,n} \mid a_{n,n+1} \end{matrix}\right. \rightarrow \left\{\begin{matrix} a'_{1,1} & a'_{1,2} & \ldots & a_{1,n} \mid a'_{1,n+1} \\ 0 & a_{2,2} & \ldots & a'_{2,n} \mid a'_{2,n+1} \\ \ldots \\ 0 & 0 & \ldots & a'_{n,n} \mid a'_{n,n+1} \end{matrix}\right.
提示

无唯一解：某行前 nn 个数均为 00。（0+0+⋯+00+0+\cdots+0）

区分无解和无数解：

  • 无解：最后结果不为 00。（0+0+⋯+0≠00+0+\cdots+0\not=0）
  • 无数解：最后结果为 00。（0+0+⋯+0=00+0+\cdots+0=0）
double a[N][N],x[N];
int sgn(double u){return (u>eps)-(u<-eps);}
int gauss(int n)
{
	int r,c;
	for(r=1,c=1;r<=n&&c<=n;r++,c++)
	{
		int t=r;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[t][c])<fabs(a[j][c]))t=j;
		}
		swap(a[r],a[t]);
		if(sgn(a[r][c])==0)
		{
			r--;
			continue;
		}
		for(int j=r+1;j<=n;j++)
		{
			if(sgn(a[j][c])==0)continue;
			double tmp=a[j][c]/a[r][c];
			for(int k=c;k<=n+1;k++)
			{
				a[j][k]-=a[r][k]*tmp;
			}
			a[j][c]=0;
		}
	} 
	for(int i=r;i<=n;i++)
	{
		if(sgn(a[i][c])!=0)return -1;
	}
	if(r<=n)return n-r+1;
	for(int i=n;i>=1;i--)
	{
		for(int j=i+1;j<=n;j++)
		{
			a[i][n+1]-=a[i][j]*x[j];
		}
		x[i]=a[i][n+1]/a[i][i];
	}
	return 0;
}

例题​

洛谷 P3389 【模板】高斯消元法​

#include <bits/stdc++.h>
using namespace std;

const double eps=1e-8;
const int N=105;
double a[N][N];
bool gauss(int n)
{
	for(int i=1;i<=n;i++)
	{
		int t=i;
		for(int j=i+1;j<=n;j++)
		{
			if(fabs(a[t][i])<fabs(a[j][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		double tmp=a[i][i];
		for(int j=1;j<=n+1;j++)a[i][j]/=tmp;
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			tmp=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;k++)
			{
				a[j][k]-=a[i][k]*tmp;
			}
		}
	}
	return 1;
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n;
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n+1;j++)
		{
			cin>>a[i][j];
		}
	}
	if(gauss(n)==0)
	{
		cout<<"No Solution"<<'\n';
		return 0;
	}
	for(int i=1;i<=n;i++)
	{
		cout<<fixed<<setprecision(2)<<a[i][n+1]<<'\n';
	}
	return 0;
}

洛谷 P2455 [SDOI2006] 线性方程组​

#include <bits/stdc++.h>
using namespace std;

const double eps=1e-8;
const int N=55;
double a[N][N],x[N];
int sgn(double x){return (x>eps)-(x<-eps);}
int gauss(int n)
{
	int r,c;
	for(r=1,c=1;r<=n&&c<=n;r++,c++)
	{
		int x=r;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[x][c])<fabs(a[j][c]))x=j;
		}
		swap(a[r],a[x]);
		if(sgn(a[r][c])==0)
		{
			r--;
			continue;
		}
		for(int j=r+1;j<=n;j++)
		{
			if(sgn(a[j][c])==0)continue;
			double tmp=a[j][c]/a[r][c];
			for(int k=c;k<=n+1;k++)
			{
				a[j][k]-=a[r][k]*tmp;
			}
			a[j][c]=0;
		}
	} 
	for(int i=r;i<=n;i++)
	{
		if(sgn(a[i][c])!=0)return -1;
	}
	if(r<=n)return n-r+1;
	for(int i=n;i>=1;i--)
	{
		for(int j=i+1;j<=n;j++)
		{
			a[i][n+1]-=a[i][j]*x[j];
		}
		x[i]=a[i][n+1]/a[i][i];
	}
	return 0;
}
int main()
{
	int n;
	cin>>n;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n+1;j++)
		{
			cin>>a[i][j];
		}
	}
	int ans=gauss(n);
	if(ans==-1)cout<<-1<<'\n';
	else if(ans>0)cout<<0<<'\n';
	else
	{
		for(int i=1;i<=n;i++)
		{
			cout<<fixed<<setprecision(2)<<'x'<<i<<'='<<x[i]<<'\n';
		}
	}
	return 0;
}

洛谷 P7112 【模板】行列式求值​

#include <bits/stdc++.h>
using namespace std;

using ll=long long;
const int inf=0x3f3f3f3f;
const int N=605;
ll a[N][N];
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n,mod;
	cin>>n>>mod;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			cin>>a[i][j];
		}
	}
	ll ans=1;
	for(int i=1;i<=n;i++)
	{
		for(int j=i+1;j<=n;j++)
		{
			while(a[i][i])
			{
	 			int tmp=a[j][i]/a[i][i];
				for(int k=i;k<=n;k++)
				{
					a[j][k]=(a[j][k]-a[i][k]*tmp%mod)%mod;
				}
				swap(a[i],a[j]);
				ans=-ans;
			}
			swap(a[i],a[j]);
			ans=-ans;
		}
	}
	for(int i=1;i<=n;i++)
	{
		ans=ans*a[i][i]%mod;
	}
	cout<<(ans+mod)%mod<<'\n';
	return 0;
}

POJ 3532 Resistance​

#include <bits/stdc++.h>
using namespace std;

const double eps=1e-8;
const int N=105;
double a[N][N];
bool gauss(int n)
{
	for(int i=1;i<=n;i++)
	{
		int t=i;
		for(int j=i+1;j<=n;j++)
		{
			if(fabs(a[t][i])<fabs(a[j][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		double tmp=a[i][i];
		for(int j=1;j<=n+1;j++)a[i][j]/=tmp;
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			tmp=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;k++)
			{
				a[j][k]-=a[i][k]*tmp;
			}
		}
	}
	return 1;
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n,m;
	cin>>n>>m;
	while(m--)
	{
		int u,v;
		double r;
		cin>>u>>v>>r;
		double g=1/r;
		a[u][u]-=g;
		a[v][v]-=g;
		a[u][v]+=g;
		a[v][u]+=g;
	}
	for(int i=1;i<=n;i++)a[n][i]=(i==1);
	a[1][n+1]=1;
	gauss(n);
	cout<<fixed<<setprecision(2)<<a[n][n+1]<<'\n';
	return 0;
}
最后于 更新