前置
- 在线性代数的课程中，我们就学过基本的矩阵及其行初等变化。
根据这些初等变化，我们的老师就高速我们怎么样去进行消元，然后求解线性方程组。
模板题
- 【模板】高斯消元法 | 洛谷 P3389
挺基础的蓝题，虽然是个模板但是思路很简单可以手敲。
⌈
\lceil
⌈高斯约旦消元法
⌋
\rfloor
⌋
- 首先，我们获得一个系数矩阵
A
n
×
n
A_{n\times n}
An×n 和一个列向量B
n
×
1
B_{n\times 1}
Bn×1
我们把它合并成一个增广矩阵A
n
×
(
n
+
1
)
A_{n\times (n+1)}
An×(n+1)
那么这个矩阵代表什么意思呢？比如我们有一下方程组：
{
1
x
+
8
y
+
2
z
=
5
2
x
+
4
y
+
12
z
=
10
1
x
+
0
y
+
4
z
=
20
\begin{cases} 1x+8y+2z=5\\ 2x+4y+12z=10\\ 1x+0y+4z=20\\ \end{cases}
⎩⎪⎨⎪⎧1x+8y+2z=52x+4y+12z=101x+0y+4z=20
我们的增广矩阵就如下所示
A
=
[
1
8
2
∣
5
2
4
12
∣
10
1
0
4
∣
20
]
A= \begin{bmatrix} 1&8&2&|&5\\ 2&4&12&|&10\\ 1&0&4&|&20\\ \end{bmatrix}
A=⎣⎡1218402124∣∣∣51020⎦⎤
【做法】：
（1）我们按列进行处理。假设现在处理到第i
i
i 列。
为了使得答案精度更高，我们选择 该列中行号≥
i
\ge i
≥i 的行中最大的值的那行 ，然后把该行和第i
i
i 行进行交换。
（2）如果某一列的所有值都为0
0
0 ，那么说明该行变量是个无关变量，任意取值都满足。
（3）根据行初等变化，我们把除了第i
i
i 行的每一行都减去一定的比值，使得 该行第i
i
i 列为0
0
0
（4）这样操作完，系数矩阵只有主对角线位置的元素非零。容易算出每一个变量的取值了。
【举例】：
A
=
[
1
8
2
∣
5
2
4
12
∣
10
1
0
4
∣
20
]
A= \begin{bmatrix} 1&8&2&|&5\\ 2&4&12&|&10\\ 1&0&4&|&20\\ \end{bmatrix}
A=⎣⎡1218402124∣∣∣51020⎦⎤
（1）我们找第一列最大的元素，为2
2
2 ，该行和第一行交换
A
=
[
2
4
12
∣
10
1
8
2
∣
5
1
0
4
∣
20
]
A= \begin{bmatrix} 2&4&12&|&10\\ 1&8&2&|&5\\ 1&0&4&|&20\\ \end{bmatrix}
A=⎣⎡2114801224∣∣∣10520⎦⎤
（2）第二行减去第一行的1
2
\frac{1}{2}
21 倍，第三行减去第一行的1
2
\frac{1}{2}
21 倍，得到新的矩阵：
A
=
[
2
4
12
∣
10
0
6
−
4
∣
0
0
−
2
−
2
∣
15
]
A= \begin{bmatrix} 2&4&12&|&10\\ 0&6&-4&|&0\\ 0&-2&-2&|&15\\ \end{bmatrix}
A=⎣⎡20046−212−4−2∣∣∣10015⎦⎤
（3）找第二列的行号为2
∼
3
2\sim 3
2∼3 的最大值。由于6
6
6 最大且本身就是第二列，因此无需交换。
（4）第一行减去第二行的4
6
\frac{4}{6}
64 倍，第三行减去第二行的−
2
6
-\frac{2}{6}
−62 倍，得到：
A
=
[
2
0
44
3
∣
10
0
6
−
4
∣
0
0
0
−
10
3
∣
15
]
A= \begin{bmatrix} 2&0&\frac{44}{3}&|&10\\ 0&6&-4&|&0\\ 0&0&-\frac{10}{3}&|&15\\ \end{bmatrix}
A=⎣⎡200060344−4−310∣∣∣10015⎦⎤
（5）找第三列的行号为3
∼
3
3\sim 3
3∼3 的最大值。最大值为−
10
3
-\frac{10}{3}
−310，无需交换。
（6）第一行减去第三行的−
44
3
⋅
3
10
-\frac{44}{3}\cdot\frac{3}{10}
−344⋅103 倍，第二行减去第三行的4
⋅
3
10
4\cdot\frac{3}{10}
4⋅103 倍。得到：
A
=
[
2
0
0
∣
76
0
6
0
∣
−
18
0
0
−
10
3
∣
15
]
A= \begin{bmatrix} 2&0&0&|&76\\ 0&6&0&|&-18\\ 0&0&-\frac{10}{3}&|&15\\ \end{bmatrix}
A=⎣⎡20006000−310∣∣∣76−1815⎦⎤
我们容易得到，变量的值为：
{
x
=
38
y
=
−
3
z
=
−
4.5
\begin{cases} x=38\\ y=-3\\ z=-4.5\\ \end{cases}
⎩⎪⎨⎪⎧x=38y=−3z=−4.5
完美！(手算错两遍咳咳咳)
补充
- 它和普通的高斯消元法有什么区别吗？
高斯消元法每次消元，最后会得到系数矩阵的上三角矩阵
然后通过回带来得到所有解，稍微麻烦些。
举个例子，比如得到下面的样子：
A
=
[
2
4
12
∣
10
0
8
2
∣
5
0
0
4
∣
20
]
A= \begin{bmatrix} 2&4&12&|&10\\ 0&8&2&|&5\\ 0&0&4&|&20\\ \end{bmatrix}
A=⎣⎡2004801224∣∣∣10520⎦⎤
然后易得z
z
z 的值，然后带回第二行得到y
y
y 的值，再带回第一行得到x
x
x 的值。
我个人感觉高斯约旦消元法更加方便。
代码
- 时间复杂度：
O
(
N
3
)
O(N^3)
O(N3)
/*
_ __ __ _ _
| | \ \ / / | | (_)
| |__ _ _ \ V /__ _ _ __ | | ___ _
| '_ \| | | | \ // _` | '_ \| | / _ \ |
| |_) | |_| | | | (_| | | | | |___| __/ |
|_.__/ \__, | \_/\__,_|_| |_\_____/\___|_|
__/ |
|___/
*/
const int MAX = 1e2+50;
const ll MOD = 1e9+7;
double aa[MAX][MAX];
bool GJ(int n){
for(int i = 1;i <= n;++i){
int now = i;
double mx = aa[i][i];
for(int j = i + 1;j <= n;++j){
if(aa[j][i] > mx){
now = j;
mx = aa[i][j];
}
}
if(fabs(mx) < EPS)return false;
if(now != i)
for(int j = 1;j <= n + 1;++j)
swap(aa[now][j],aa[i][j]);
for(int j = 1;j <= n;++j){
if(j == i)continue;
for(int k = n + 1;k >= i;--k){
aa[j][k] -= aa[j][i] / aa[i][i] * aa[i][k];
}
}
}
return true;
}
int main()
{
int n;cin >> n;
for(int i = 1;i <= n;++i)
for(int j = 1;j <= n + 1;++j)
scanf("%lf",&aa[i][j]);
bool can = GJ(n);
if(!can)puts("No Solution");
else{
for(int i = 1;i <= n;++i)
printf("%.2f\n",aa[i][n+1] / aa[i][i]);
}
return 0;
}
更新
- 突然出了这么一题 线性方程组 | 洛谷 P2455
除了输出唯一解外，还需要判定无解与无穷解的情况。
虽然高斯约旦消元法不大擅长判断这种无解还是无穷解的情况，但是我们仍然可以用高斯约旦消元法去做。
主要思路就是我们设一个当前行r
r
r，如果你找不到当前列的主元(就是全是0)，那么直接跳到下一列(但是当前行不变)；否则，就正常的消元，然后当前行递增。
当你消元完毕后，如果当前行r
<
n
r<n
r<n，那么表示最后几行的系数矩阵全为0
0
0。我们需要判断这几行的列向量是否全0
0
0，以此判断是无解还是无穷解。
改进代码
double aa[MAX][MAX]; /// 增广矩阵
double ans[MAX];
int GJ(int n){
int r = 1;
for(int c = 1;c <= n;++c){
int now = r;
double mx = fabs(aa[r][c]); /// 找最大系数改为找绝对值最大系数
for(int i = r + 1;i <= n;++i){
if(fabs(aa[i][c]) > mx){
now = i;
mx = fabs(aa[i][c]);
}
}
if(mx < EPS)continue; /// 跳过当前列
for(int j = 1;j <= n + 1;++j)
swap(aa[r][j],aa[now][j]);
for(int i = 1;i <= n;++i){
if(i == r)continue;
for(int j = c + 1;j <= n + 1;++j)
aa[i][j] -= aa[i][c] / aa[r][c] * aa[r][j];
}
r++; /// 当前行变为下一行
}
r--;
if(r < n){
for(int i = r + 1;i <= n;++i)
if(fabs(aa[i][n + 1]) > EPS)return -1; /// 无解
return 0; /// 无穷解
}
for(int i = 1;i <= n;++i){
ans[i] = aa[i][n + 1] / aa[i][i];
if(fabs(ans[i]) < EPS)ans[i] = 0; /// 注意把输出 -0改成 0
}
return 1; /// 有唯一解
}