「学习笔记」高斯消元
用于求解线性方程组
「学习笔记」高斯消元
点击查看目录
注意：本文内的高斯消元均是「高斯-约旦消元法（\(\text{Gauss-Jordan Elimination}\)）」。
突然发现有的电脑好像看不了增广矩阵里的竖线，那就凑合看吧（
算法
思路
我们举个例子：
\[\begin{cases}
x_2+x_3&=3\\
x_1+2x_3&=9\\
2x_1+x_2+x_3&=13
\end{cases}
\]
我们把它的所有系数化为矩阵 \(A\)：
\[\left[\begin{matrix}
0&1&1\\
1&0&2\\
2&1&1\\
\end{matrix}
\middle|
\begin{matrix}
3 \\
9 \\
13
\end{matrix}
\right]
\]
然后开始推：
- 首先遍历第 \(i\) 列的第 \(i\) 行至 \(n\) 行（因为前 \(i-1\) 行已经形成阵型了，不能破坏），找出这一列最大的系数，并将那一行交换至第 \(i\) 行。（当前 \(i=1\)）
\[\left[\begin{matrix}
2&1&1\\
1&0&2\\
0&1&1\\
\end{matrix}
\middle|
\begin{matrix}
13 \\
9 \\
3 \\
\end{matrix}
\right]
\]
-
然后判断解，即判断 \(A_{i,i}\) 是否为零。若为零则无唯一解，否则有唯一解。（这个例子显然有唯一解）
-
系数化为一，即将 \(A_{i,j}(i<j\le n+1)\) 除以 \(A_{i,i}\)，再将 \(A_{i,i}\) 除以自己得 \(1\)。
\[\left[\begin{matrix}
1&0.5&0.5\\
1&0&2\\
0&1&1\\
\end{matrix}
\middle|
\begin{matrix}
6.5 \\
9 \\
3 \\
\end{matrix}
\right]
\]
- 然后用加减消元法消掉一部分元：
\[\left[\begin{matrix}
1&0.5&0.5\\
0&-0.5&1.5\\
0&1&1\\
\end{matrix}
\middle|
\begin{matrix}
6.5 \\
2.5 \\
3 \\
\end{matrix}
\right]
\]
- 重复上述步骤：
\[\left[\begin{matrix}
1&0.5&0.5\\
0&-0.5&1.5\\
0&1&1\\
\end{matrix}\middle|\begin{matrix}
6.5 \\
2.5 \\
3 \\
\end{matrix}\right]
\Longrightarrow
\left[\begin{matrix}
1&0.5&0.5\\
0&1&1\\
0&-0.5&1.5\\
\end{matrix}\middle|\begin{matrix}
6.5 \\
3 \\
2.5 \\
\end{matrix}\right]
\Longrightarrow
\left[\begin{matrix}
1&0&0\\
0&1&1\\
0&0&2\\
\end{matrix}\middle|\begin{matrix}
5 \\
3 \\
4 \\
\end{matrix}\right]
\Longrightarrow
\left[\begin{matrix}
1&0&0\\
0&1&1\\
0&0&1\\
\end{matrix}\middle|\begin{matrix}
5 \\
3 \\
2 \\
\end{matrix}\right]
\Longrightarrow
\left[\begin{matrix}
1&0&0\\
0&1&0\\
0&0&1\\
\end{matrix}\middle|\begin{matrix}
5 \\
1 \\
2 \\
\end{matrix}\right]
\]
所有操作步骤完成，最后得到解：
\[\begin{cases}
x_1=5\\
x_2=1\\
x_3=2\\
\end{cases}
\]
时间复杂度 \(\Theta(n^3)\)。
注意：输入不要用快读，虽然有的题输入只有整数，但很多题都不是！
代码
点击查看代码
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
_for(j,1,n+1)swap(a[i][j],a[mx][j]);
if(fabs(a[i][i])<eps){puts("No Solution");return;}
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1;
_for(j,1,n){
if(i==j)continue;
_for(k,i+1,n+1)a[j][k]-=a[j][i]*a[i][k];
}
}
_for(i,1,n)printf("%.2lf\n",a[i][n+1]);
}
例题:[SDOI2006]线性方程组
思路
本题主要难点在于如何区分无解和无限组解。
举个例子：
\[A=
\left[\begin{matrix}
1&0&1\\
0&2&1\\
1&2&2\\
\end{matrix}\middle|\begin{matrix}
2 \\
3 \\
5 \\
\end{matrix}\right],
B=
\left[\begin{matrix}
1&0&1\\
0&2&1\\
1&2&2\\
\end{matrix}\middle|\begin{matrix}
2 \\
3 \\
7 \\
\end{matrix}\right],
C=
\left[\begin{matrix}
1&0&1\\
0&2&1\\
1&1&2\\
\end{matrix}\middle|\begin{matrix}
2 \\
3 \\
5 \\
\end{matrix}\right]
\]
-
在 \(A\) 中，第三行的前 \(3\) 项是前两行相加的结果，第 \(4\) 项也是，所以有用的行只有两个（这个数量称为秩），有一个数是不确定的，所以有不唯一解。
-
在 \(B\) 中，第三行的前 \(3\) 项是前两行相加的结果，但第 \(4\) 项不是，会出现 \(0=2\) 的情况，所以无解。
-
在 \(C\) 中，它的秩是 \(3\)，所以有唯一解：
\[\begin{cases}
x_1=-1\\
x_2=0\\
x_3=3\\
\end{cases}
\]
归纳一下：
-
若秩不等于未知数的数量且消元后出现 \(\sum_{j=1}^{n}a_{i,j}=0\) 且 \(a_{i,n+1}\neq0\) 的情况，则无解。
-
若秩不等于未知数的数量且消元后出现 \(\sum_{j=1}^{n+1}a_{i,j}=0\) 的情况，则有无数组解。
-
否则有唯一解。
代码
点击查看代码
const ll N=110,inf=1ll<<40;
const double eps=1e-9;
ll n,l=1;db a[N][N];
namespace SOLVE{
inline void Gauss(){
_for(i,1,n){
ll mx=l;
_for(j,l+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
if(fabs(a[mx][i])<eps)continue;
swap(a[mx],a[l]);
_for(j,i+1,n+1)a[l][j]/=a[l][i];
a[l][i]=1;
_for(j,1,n){
if(j==l)continue;
_for(k,i+1,n+1)a[j][k]-=a[l][k]*a[j][i];
}
++l;
}
}
inline void Jie(){
if(l<=n){
_for(i,l,n)if(fabs(a[i][n+1])>eps){puts("-1");return;}
puts("0");return;
}
_for(i,1,n){
printf("x%lld=%.2lf\n",i,a[i][n+1]);
}
}
inline void In(){
scanf("%lld",&n);
_for(i,1,n)_for(j,1,n+1)scanf("%lf",&a[i][j]);
Gauss();
Jie();
return;
}
}
练习题
[JSOI2008]球形空间产生器
思路
设球心为 \(p(p_1,p_2,\cdots,p_n)\)，可列式子：
\[\begin{aligned}
r=&(a_{1,1}-p_1)^2+(a_{1,2}-p_2)^2+\cdots+(a_{1,n}-p_n)^2\\
=&(a_{2,1}-p_1)^2+(a_{2,2}-p_2)^2+\cdots+(a_{2,n}-p_n)^2\\
=&\cdots\\
=&(a_{n+1,1}-p_1)^2+(a_{n+1,2}-p_2)^2+\cdots+(a_{n+1,n}-p_n)^2\\
\end{aligned}
\]
拆一下得：
\[\begin{aligned}
r=&(a_{1,1}^2-2a_{1,1}p_1+p_1^2)+(a_{1,2}^2-2a_{1,2}p_2+p_2^2)+\cdots+(a_{1,n}^2-2a_{1,n}p_n+p_n^2)\\
=&(a_{2,1}^2-2a_{2,1}p_1+p_1^2)+(a_{2,2}^2-2a_{2,2}p_2+p_2^2)+\cdots+(a_{2,n}^2-2a_{2,n}p_n+p_n^2)\\
=&\cdots\\
=&(a_{n+1,1}^2-2a_{n+1,1}p_1+p_1^2)+(a_{n+1,2}^2-2a_{n+1,2}p_2+p_2^2)+\cdots+(a_{n+1,n}^2-2a_{n+1,n}p_n+p_n^2)\\
\end{aligned}
\]
化简得：
\[\begin{aligned}
r=&\sum_{i=1}^{n}a_{1,i}^2-\sum_{i=1}^{n}2a_{1,i}p_i\\
=&\sum_{i=1}^{n}a_{2,i}^2-\sum_{i=1}^{n}2a_{2,i}p_i\\
=&\cdots\\
=&\sum_{i=1}^{n}a_{n+1,i}^2-\sum_{i=1}^{n}2a_{n+1,i}p_i\\
\end{aligned}
\]
即：
\[\begin{aligned}
\sum_{i=1}^{n}a_{1,i}^2=&r+\sum_{i=1}^{n}2a_{1,i}p_i\\
\sum_{i=1}^{n}a_{2,i}^2=&r+\sum_{i=1}^{n}2a_{2,i}p_i\\
\cdots=&r+\cdots\\
\sum_{i=1}^{n}a_{n+1,i}^2=&r+\sum_{i=1}^{n}2a_{n+1,i}p_i\\
\end{aligned}
\]
发现就是个方程，但是这个 \(r\) 混在里面很难受。那么我们把它当成一个系数为 \(1\) 的未知数，则原式可以化成如下增广矩阵：
\[\left[
\begin{matrix}
2a_{1,1}&2a_{1,2}&\cdots&2a_{1,n}&1\\\\
2a_{2,1}&2a_{2,2}&\cdots&2a_{2,n}&1\\\\
\cdots&\cdots&\cdots&\cdots&1\\\\
2a_{n+1,1}&2a_{n+1,2}&\cdots&2a_{n+1,n}&1
\end{matrix}
\middle|
\begin{matrix}
\sum_{i=1}^{n}a_{1,i}^2\\\\
\sum_{i=1}^{n}a_{2,i}^2\\\\
\cdots\\\\
\sum_{i=1}^{n}a_{n+1,i}^2\\
\end{matrix}
\right]
\]
就可以直接高斯消元了。
Latex 当算草纸真好用。
代码
点击查看代码
const ll N=20,inf=1ll<<40;
const double eps=1e-6;
ll n;ldb a[N][N];
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
swap(a[mx],a[i]);
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1.0;
_for(j,1,n){
if(i==j)continue;
for_(k,n+1,i)a[j][k]-=a[j][i]*a[i][k];
}
}
}
inline void In(){
scanf("%lld",&n),++n;
_for(i,1,n){
_for(j,1,n-1){
scanf("%Lf",&a[i][j]);
a[i][n+1]+=a[i][j]*a[i][j];
a[i][j]*=2.0;
}
a[i][n]=1;
}
Gauss();
_for(i,1,n-1)printf("%.3Lf ",a[i][n+1]);
return;
}
}
[USACO10HOL]Driving Out the Piggies G
思路
好题，还对高斯消元能出什么题没什么认知就 刷新了我对高斯消元能出什么题的认知。
我们设到节点 \(i\) 的期望步数为 \(x_i\)，节点 \(i\) 在图上的出度为 \(d_i\)。
那么显然：
\[x_i=[i=1]+\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times \frac{1}{d_j}
\]
显然 \(x_i\times\frac{p}{q}\) 为实际爆炸的概率，那么考虑如何求 \(x_i\)。我们进行一个移项：
\[x_i-\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times \frac{1}{d_j}=[i=1]
\]
直接冲一个高斯消元即可。
代码
点击查看代码
const ll N=310,inf=1ll<<40;
const double eps=1e-13;
ll n,m,tu[N][N];ldb p,q,d[N],a[N][N];
namespace SOLVE{
inline ll rnt(){
ll x=0,w=1;char c=getchar();
while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x*w;
}
inline void Pre(){
a[1][n+1]=1;
_for(i,1,n){
_for(j,1,n){
if(i==j)a[i][j]=1.0;
else if(tu[i][j])a[i][j]=-(1.0-p/q)/d[j];
}
}
}
inline void Gauss(){
_for(i,1,n){
ll mx=i;
_for(j,i+1,n)if(a[mx][i]<a[j][i])mx=j;
swap(a[mx],a[i]);
_for(j,i+1,n+1)a[i][j]/=a[i][i];
a[i][i]=1;
_for(j,1,n){
if(i==j)continue;
for_(k,n+1,i)a[j][k]-=a[i][k]*a[j][i];
}
}
}
inline void In(){
scanf("%lld%lld%Lf%Lf",&n,&m,&p,&q);
_for(i,1,m){
ll x=rnt(),y=rnt();
d[x]+=1.0,d[y]+=1.0;
tu[x][y]=tu[y][x]=1;
}
Pre(),Gauss();
_for(i,1,n)printf("%.9Lf\n",a[i][n+1]*p/q);
return;
}
}
\[\Huge{\mathfrak{The\ End}}
\]
