洛谷端文章：https://www.luogu.com.cn/article/n5rjrsdw，如无法访问访问https://www.luogu.me/article/n5rjrsdw。
理论
这个算法可以求 n n n 元一次方程组的解。
举个例子，现在已知有三元一次方程：
{ x + 2 y + 3 z = 36 − x + y + 2 z = 17 x + 2 z = 17 \begin{cases} x + 2y + 3z = 36 \\ -x + y + 2z = 17 \\ x + 2z = 17 \end{cases} ⎩⎪⎨⎪⎧x+2y+3z=36−x+y+2z=17x+2z=17
解得：
{
x
=
3
y
=
6
z
=
7
\begin{cases} x = 3 \\ y = 6 \\ z = 7 \end{cases}
⎩⎪⎨⎪⎧x=3y=6z=7
将问题扩展为
n
n
n 个未知数的方程组：
{
a
1
,
1
x
1
+
a
1
,
2
x
2
+
⋯
+
a
1
,
n
x
n
=
b
1
a
2
,
1
x
1
+
a
2
,
2
x
2
+
⋯
+
a
2
,
n
x
n
=
b
2
⋮
a
n
,
1
x
1
+
a
n
,
2
x
2
+
⋯
+
a
n
,
n
x
n
=
b
n
\begin{cases} a_{1,1}x_1 + a_{1,2}x_2 + \cdots + a_{1,n}x_n = b_1 \\ a_{2,1}x_1 + a_{2,2}x_2 + \cdots + a_{2,n}x_n = b_2 \\ \vdots \\ a_{n,1}x_1 + a_{n,2}x_2 + \cdots + a_{n,n}x_n = b_n \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧a1,1x1+a1,2x2+⋯+a1,nxn=b1a2,1x1+a2,2x2+⋯+a2,nxn=b2⋮an,1x1+an,2x2+⋯+an,nxn=bn
对应的矩阵表示为：
[
a
1
,
1
a
1
,
2
⋯
a
1
,
n
b
1
a
2
,
1
a
2
,
2
⋯
a
2
,
n
b
2
⋮
⋮
⋱
⋮
⋮
a
n
,
1
a
n
,
2
⋯
a
n
,
n
b
n
]
\begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} & b_1 \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} & b_n \\ \end{bmatrix}
⎣⎢⎢⎢⎡a1,1a2,1⋮an,1a1,2a2,2⋮an,2⋯⋯⋱⋯a1,na2,n⋮an,nb1b2⋮bn⎦⎥⎥⎥⎤
高斯消元步骤分解
我们以刚才的例子。
初始矩阵：
[ 1 2 3 36 − 1 1 2 17 1 0 2 17 ] \begin{bmatrix} 1 & 2 & 3 & 36 \\ -1 & 1 & 2 & 17 \\ 1 & 0 & 2 & 17 \\ \end{bmatrix} ⎣⎡1−11210322361717⎦⎤
- 处理第一列
消去第二、三行的第一列：
[ 1 2 3 36 0 3 5 53 0 − 2 − 1 − 19 ] \begin{bmatrix} 1 & 2 & 3 & 36 \\ 0 & 3 & 5 & 53 \\ 0 & -2 & -1 & -19 \\ \end{bmatrix} ⎣⎡10023−235−13653−19⎦⎤
- 处理第二列
第二行归一化并消去其他行：
[ 1 0 − 1 3 14 3 0 1 5 3 53 3 0 0 7 3 49 3 ] \begin{bmatrix} 1 & 0 & -\frac{1}{3} & \frac{14}{3} \\ 0 & 1 & \frac{5}{3} & \frac{53}{3} \\ 0 & 0 & \frac{7}{3} & \frac{49}{3} \\ \end{bmatrix} ⎣⎡100010−313537314353349⎦⎤
- 处理第三列：
第三行归一化并回代：
[ 1 0 0 3 0 1 0 6 0 0 1 7 ] \begin{bmatrix} 1 & 0 & 0 & 3 \\ 0 & 1 & 0 & 6 \\ 0 & 0 & 1 & 7 \\ \end{bmatrix} ⎣⎡100010001367⎦⎤
具体操作步骤
- 对每个 x i x_i xi，选取系数绝对值最大的一行，交换至第 i i i 行并归一化。
- 用主行消去其他行中 i i i 的系数。
- 回代：从最后一行倒序求解变量值，公式为：
时间复杂度 O ( n 3 ) \text{O}(n^3) O(n3)。
特殊情况分析
无解情况
形如这样的，方程一样但是答案不一样，比如：
[ 3 7 5 10 3 7 5 20 ] \begin{bmatrix} 3 & 7 & 5 & 10 \\ 3 & 7 & 5 & 20 \\ \end{bmatrix} [3377551020]
是不行的。
还有形如这样的：
[
0
0
0
1
−
4
4
1
13
]
\begin{bmatrix} 0 & 0 & 0 & 1 \\ -4 & 4 & 1 & 13 \\ \end{bmatrix}
[0−40401113]
存在矛盾方程
0
x
+
0
y
+
0
z
≠
0
0x+0y+0z \ne 0
0x+0y+0z=0，方程组无解。
一般情况：一整行的系数都为 0 0 0，答案不为 0 0 0。
无穷解情况
[
0
1
2
18
0
6
−
3
4
]
\begin{bmatrix} 0 & 1 & 2 & 18 \\ 0 & 6 & -3 & 4 \\ \end{bmatrix}
[00162−3184]
第一列全零，自由变量存在，方程组有无穷多解。
形如这样的，一整列的系数都为 0 0 0，说明这个数没有参与任何方程。是多少都无所谓，与答案没有关联。
模板
P2455 [SDOI2006] 线性方程组
#include<bits/stdc++.h>
using namespace std;
const int N=160;
const double eps=1e-6;
double a[N][N];
int n;
int gauss()
{
int c,r;//C:当前处理X_c,R:处理了 r 个未知数
for(c=0,r=0;c<n;c++)
{
int t=r;
for(int i=r;i<n;i++)
{
if(fabs(a[i][c])>fabs(a[t][c]))//如果这个系数更大
{
t=i;//更新
}
}
if(fabs(a[t][c])<eps)
{//如果最大的数还是 0，则代表所有方程这个项的系数都是 0。
continue;//就不管
}
for(int i=c;i<n+1;i++)
{
swap(a[r][i],a[t][i]);//将其交换至第 c 行
}
for(int i=n;i>=c;i--)
{
a[r][i]/=a[r][c];//除以系数
}
for(int i=0;i<n;i++)//消消乐
{
if(fabs(a[i][c])>eps && i!=r)//如果这一项系数不为 0 且不是系数最大的一行
{
for(int j=n;j>=c;j--)
{
a[i][j]-=a[i][c]*a[r][j];//消元
}
}
}
r++;//新处理了一个未知数
}
if(r<n)//如果没有处理完所有的未知数
{
for(int i=r;i<n;i++)//检查剩下的方程
{
if(fabs(a[i][n])>eps)
{//如果有最终答案>=0说明无解
return -1;
}
}
return 0;//反之无解
}
return 1;//有唯一解
}
int main()
{
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n;
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
{
cin>>a[i][j];
}
}
int ans=gauss();
if(ans==1)
{
for(int i=0;i<n;i++)
{
printf("x%d=%.2f\n",i+1,a[i][n]);
}
}
else
{
cout<<ans;
}
return 0;
}
P3389 【模板】高斯消元法
这个
gauss
\operatorname{gauss}
gauss 是经典void类型的。
#include<bits/stdc++.h>
using namespace std;
const int MAXN=100;
int n;
double a[MAXN+5][MAXN+5];
void gauss(int x)
{
for(int i=0;i<x;i++)
{
int row=i;
for(int j=i;j<x;j++)
{
if(fabs(a[row][i])<fabs(a[j][i]))
{
row=j;
}
}
if(row!=i)
{
swap(a[row],a[i]);
}
double div1=a[i][i];
for(int j=0;j<=x;j++)
{
a[i][j]/=div1;
}
for(int j=i+1;j<x;j++)
{
double div2=a[j][i];
for(int k=0;k<=x;k++)
{
a[j][k]-=a[i][k]*div2;
}
}
}
//-----------------------------复原
for(int i=x-1;i>=0;i--)
{
for(int j=i-1;j>=0;j--)
{
a[j][x]-=a[j][i]*a[i][x];
a[j][i]=0;
}
}
//------------------------------
}
int main()
{
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n;
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
{
cin>>a[i][j];
}
}
gauss(n);
for(int i=0;i<n;i++)
{
bool flag=1;
for(int j=0;j<n;j++)
{
if(fabs(a[i][j])>1e-8)
{
flag=0;
}
}
if(flag)
{
cout<<"No Solution";
return 0;
}
}
for(int i=0;i<n;i++)
{
cout<<fixed<<setprecision(2)<<a[i][n]<<'\n';
}
return 0;
}
扫一扫
405
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
