Skip to main content

消元

参考资料​

简介​

高斯消元法（Gaussian Elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。

{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+\dots+a_{1,n}x_n=a_{1,n+1} \\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n=a_{2,n+1} \\ \vdots \\ a_{n,1}x_1+a_{n,2}x_2+\dots+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} & \dots & a_{1,n} & \mid a_{1,n+1} \\ a_{2,1} & a_{2,2} & \dots & a_{2,n} & \mid a_{2,n+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} & \mid a_{n,n+1} \end{matrix} \right. \to \left\{ \begin{matrix} 1 & 0 & \dots & 0 & \mid a'_{1,n+1} \\ 0 & 1 & \dots & 0 & \mid a'_{2,n+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & \mid a'_{n,n+1} \end{matrix} \right. 
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[j][i])>fabs(a[t][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		for(int j=n+1;j>=i;j--)a[i][j]/=a[i][i];
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			for(int k=n+1;k>i;k--)a[j][k]-=a[i][k]*a[j][i];
		}
	}
	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} & \dots & a_{1,n} & \mid a_{1,n+1} \\ a_{2,1} & a_{2,2} & \dots & a_{2,n} & \mid a_{2,n+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} & \mid a_{n,n+1} \end{matrix} \right. \to \left\{ \begin{matrix} a'_{1,1} & a'_{1,2} & \dots & a_{1,n} & \mid a'_{1,n+1} \\ 0 & a_{2,2} & \dots & a'_{2,n} & \mid a'_{2,n+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & a'_{n,n} & \mid a'_{n,n+1} \end{matrix} \right.

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

  • 无解：最后结果不为 00。（0+0+⋯+0≠00+0+\dots+0\ne 0）
  • 无穷多解：最后结果为 00。（0+0+⋯+0=00+0+\dots+0=0）
int gauss(int n)
{
	int r,c;
	for(r=1,c=1;r<=n&&c<=n;c++)
	{
		int t=r;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[j][c])>fabs(a[t][c]))t=j;
		}
		swap(a[r],a[t]);
		if(fabs(a[r][c])<eps)continue;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[j][c])<eps)continue;
			double f=a[j][c]/a[r][c];
			for(int k=c;k<=n+1;k++)a[j][k]-=a[r][k]*f;
			a[j][c]=0;
		}
		r++;
	}
	for(int i=r;i<=n;i++)
	{
		if(fabs(a[i][n+1])>eps)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]*a[j][0];
		a[i][0]=a[i][n+1]/a[i][i];
	}
	return 0;
}

行列式求值​

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

using ll=long long;
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 f=a[j][i]/a[i][i];
				for(int k=i;k<=n;k++)
				{
					a[j][k]=(a[j][k]-a[i][k]*f%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;
}

例题​

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

给定一个线性方程组，对其求解。

{a1,1x1+a1,2x2+⋯+a1,nxn=b1a2,1x1+a2,2x2+⋯+a2,nxn=b2…an,1x1+an,2x2+⋯+an,nxn=bn\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \dots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \dots + a_{2, n} x_n = b_2 \\ \dots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \dots + a_{n, n} x_n = b_n \end{cases}
Code (1)
#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[j][i])>fabs(a[t][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		for(int j=n+1;j>=i;j--)a[i][j]/=a[i][i];
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			for(int k=n+1;k>i;k--)a[j][k]-=a[i][k]*a[j][i];
		}
	}
	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] 线性方程组

已知 nn 元线性一次方程组。

{a1,1x1+a1,2x2+⋯+a1,nxn=b1a2,1x1+a2,2x2+⋯+a2,nxn=b2…an,1x1+an,2x2+⋯+an,nxn=bn\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \dots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \dots + a_{2, n} x_n = b_2 \\ \dots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \dots + a_{n, n} x_n = b_n \end{cases}

请根据输入的数据，编程输出方程组的解的情况。

Code (1)
#include <bits/stdc++.h>
using namespace std;

const double eps=1e-8;
const int N=55;
double a[N][N];
int gauss(int n)
{
	int r,c;
	for(r=1,c=1;r<=n&&c<=n;c++)
	{
		int t=r;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[j][c])>fabs(a[t][c]))t=j;
		}
		swap(a[r],a[t]);
		if(fabs(a[r][c])<eps)continue;
		for(int j=r+1;j<=n;j++)
		{
			if(fabs(a[j][c])<eps)continue;
			double f=a[j][c]/a[r][c];
			for(int k=c;k<=n+1;k++)a[j][k]-=a[r][k]*f;
			a[j][c]=0;
		}
		r++;
	}
	for(int i=r;i<=n;i++)
	{
		if(fabs(a[i][n+1])>eps)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]*a[j][0];
		a[i][0]=a[i][n+1]/a[i][i];
	}
	return 0;
}
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];
		}
	}
	int ans=gauss(n);
	if(ans==-1){cout<<-1<<'\n';return 0;}
	if(ans>0){cout<<0<<'\n';return 0;}
	for(int i=1;i<=n;i++)
	{
		cout<<fixed<<setprecision(2)<<'x'<<i<<'='<<a[i][0]<<'\n';
	}
	return 0;
}
洛谷 P7112 【模板】行列式求值

给定一个 nn 阶行列式 AA，求 ∣A∣|A|。结果对 pp 取模。

Code (1)
#include <bits/stdc++.h>
using namespace std;

using ll=long long;
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 f=a[j][i]/a[i][i];
				for(int k=i;k<=n;k++)
				{
					a[j][k]=(a[j][k]-a[i][k]*f%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;
}
SP2415 RESIST - Kirchhof Law

给定一个有 nn 个点，mm 条边的无向图（保证连通），每条边上各有一个阻值不同的电阻，求点 11 到点 nn 的总电阻，n≤100,m≤300n\le 100,m\le 300，注意是多测。

Solution

题意简述​

给定一张 nn 个节点 mm 条边的无向图，边权表示电阻阻值，求节点 11 到 nn 的等效电阻。

参考资料​

解题思路​

设节点 xx 的电势为 φx\varphi_x，边 (u,v)(u,v) 的电阻为 Ru,vR_{u,v}，则电导为：

Gu,v=1Ru,vG_{u,v}=\frac{1}{R_{u,v}}

根据 欧姆定律，从 uu 到 vv 的电流为：

Iu,v=φu−φvRu,v=Gu,v(φu−φv)I_{u,v}=\frac{\varphi_u-\varphi_v}{R_{u,v}}=G_{u,v}(\varphi_u-\varphi_v)

根据 基尔霍夫电流定律，对于任意一个普通节点，所有关联支路电流的代数和为 00。

考虑节点 uu 的电流：

∑(u,v)∈EIu,v=∑(u,v)∈EGu,v(φu−φv)=0\sum_{(u,v)\in E}I_{u,v}=\sum_{(u,v)\in E}G_{u,v}(\varphi_u-\varphi_v)=0

展开可得：

(∑(u,v)∈EGu,v)φu−∑(u,v)∈EGu,vφv=0\left(\sum_{(u,v)\in E}G_{u,v}\right)\varphi_u-\sum_{(u,v)\in E}G_{u,v}\varphi_v=0

这是一个关于各点电势的线性方程。

au,1φ1+au,2φ2+⋯+au,nφn=bua_{u,1}\varphi_1+a_{u,2}\varphi_2+\dots+a_{u,n}\varphi_n=b_u

考虑添加一条边 (u,v)(u,v)，设其电导为 GG。节点 uu 的方程加入 G(φu−φv)G(\varphi_u-\varphi_v)，节点 vv 的方程加入 G(φv−φu)G(\varphi_v-\varphi_u)，因此对应到系数矩阵就是：

au,u←au,u+Ga_{u,u}\gets a_{u,u}+G au,v←au,v−Ga_{u,v}\gets a_{u,v}-G av,v←av,v+Ga_{v,v}\gets a_{v,v}+G av,u←av,u−Ga_{v,u}\gets a_{v,u}-G

把所有边的贡献都加到矩阵里，就得到了整个线性方程组。

我们可以在节点 11 通入 1A1\mathrm{A} 电流，并将节点 nn 设为零电势点，求出节点 11 的电势。

b1=1,φn=0b_1=1,\varphi_n=0

利用 高斯消元 求解各点电势后，等效电阻为：

R=UI=φ1−φn1=φ1R=\frac{U}{I}=\frac{\varphi_1-\varphi_n}{1}=\varphi_1

参考代码​

#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[j][i])>fabs(a[t][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		for(int j=n+1;j>=i;j--)a[i][j]/=a[i][i];
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			for(int k=n+1;k>i;k--)a[j][k]-=a[i][k]*a[j][i];
		}
	}
	return 1;
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n,m;
	while(cin>>n>>m)
	{
		while(m--)
		{
			int u,v;
			double r;
			cin>>u>>v>>r;
			double g=1/r;
			a[u][u]+=g;
			a[u][v]-=g;
			a[v][v]+=g;
			a[v][u]-=g;
		}
		a[1][n+1]=1;
		for(int i=1;i<n;i++)a[n][i]=0;
		a[n][n]=1;
		a[n][n+1]=0;
		gauss(n);
		cout<<fixed<<setprecision(2)<<a[1][n+1]<<'\n';
		memset(a,0,sizeof a);
	}
	return 0;
}
POJ 3532 Resistance

给定一张图，表示一个电阻电路，求等效电阻。

Code (1)
#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[j][i])>fabs(a[t][i]))t=j;
		}
		if(fabs(a[t][i])<eps)return 0;
		swap(a[i],a[t]);
		for(int j=n+1;j>=i;j--)a[i][j]/=a[i][i];
		for(int j=1;j<=n;j++)
		{
			if(i==j)continue;
			for(int k=n+1;k>i;k--)a[j][k]-=a[i][k]*a[j][i];
		}
	}
	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[u][v]-=g;
		a[v][v]+=g;
		a[v][u]-=g;
	}
	a[1][n+1]=1;
	for(int i=1;i<n;i++)a[n][i]=0;
	a[n][n]=1;
	a[n][n+1]=0;
	gauss(n);
	cout<<fixed<<setprecision(2)<<a[1][n+1]<<'\n';
	return 0;
}