高斯消元
题目1
洛谷 P3389 【模板】高斯消元法
总的来说，就是求解一个 \(n\) 元一次方程组。
高斯消元
思路：
首先把所有系数看成一个矩阵：
\[\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5} \\
\\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5} \\
\\
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} \\
\\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} \\
\\
a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} \\
\end{bmatrix}
\]
然后将所有方程等号右边的常数放到最后一列：
\[\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & a_{1,5} & b_1 \\
\\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & a_{2,5} & b_2 \\
\\
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} & b_3 \\
\\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} & b_4 \\
\\
a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} & b_5 \\
\end{bmatrix}
\]
这样就构成了增广矩阵。
步骤：
- 枚举每一行 \(i\)，并将这一行的第 \(i\) 个未知数作为这一行的主元 \(x_i\)。 如果这个未知数的系数为 \(0\)，那就从这行开始，向下几行寻找一个该未知数系数不为 \(0\) 的方程，找到了就将这一行与当前行交换（找不到的情况就是 \(No \space Solution\)，为什么后面会说）。
- 下来再将该主元的系数化为 \(1\)，用这个主元去消下面几个方程的相同未知数的系数，即：把矩阵中与该主元在同一列且在该主元下面的系数消为 \(0\)。
- 让后我们就可以从下往上带入求解。
完成 \(1、2\) 步后，我们就会得到一个除最后一个“常数列”外的上三角矩阵，即：除了最后一列，以这个矩阵的 “左上—右下” 对角线为界，其左下角（不包括对角线）都为 \(0\)。
如下图：
\[\begin{bmatrix}
a_{1,1}' & a_{1,2}' & a_{1,3}' & a_{1,4}' & a_{1,5}' & b_1' \\
\\
0 & a_{2,2}' & a_{2,3}' & a_{2,4}' & a_{2,5}' & b_2' \\
\\
0 & 0 & a_{3,3}' & a_{3,4}' & a_{3,5}' & b_3' \\
\\
0 & 0 & 0 & a_{4,4}' & a_{4,5}' & b_4' \\
\\
0 & 0 & 0 & 0 & a_{5,5}' & b_5' \\
\end{bmatrix}
\]
最后，从下往上带入求解每一元。
代码
#include <bits/stdc++.h>
#define mkpr make_pair
#define fir first
#define sec second
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
const int maxn = 1e2 + 7;
const double eps = 1e-9;
int n;
double a[maxn][maxn];
bool gauss() {
for (int i = 1; i <= n; ++i) {
int r = i;
for (int j = i; j <= n; ++j)
if (abs(a[r][i]) < abs(a[j][i])) // 这里找该元系数最大的方程可以减小误差
r = j;
if (r != i) swap(a[r], a[i]);
if (abs(a[i][i]) < eps) return 0;
for (int j = n + 1; j >= i; --j)
a[i][j] /= a[i][i];
for (int j = i + 1; j <= n; ++j)
for (int k = n + 1; k >= i; --k)
a[j][k] -= a[i][k] * a[j][i];
}
for (int i = n - 1; i >= 1; --i)
for (int j = i + 1; j <= n; ++j)
a[i][n + 1] -= a[i][j] * a[j][n + 1];
return 1;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n + 1; ++j)
scanf("%lf", &a[i][j]);
if (gauss()) {
for (int i = 1; i <= n; ++i)
printf("%.2lf\n", a[i][n + 1]);
} else {
printf("No Solution\n");
}
return 0;
}
高斯-约旦消元
思路
同上，将系数、常数视为一个矩阵。不同的是如何消元。
步骤：
- 依旧是枚举 \(i\) 行，循环找到主元 \(x_i\) 的系数不为 \(0\) 的方程。
- 再次拿 \(x_i\) 去消其他所有相同的元，“消去所有” 指的就是矩阵中与 \(a_{i,i}\)（即此方程中 \(x_i\) 的系数） 在同一列的、除 \(a_{i,i}\) 自己的其他系数化为 \(0\)。
最后就得到了一条只有对角线的系数矩阵：
\[\begin{bmatrix}
a_{1,1}' & 0 &0 & 0 & 0 & b_1' \\
\\
0 & a_{2,2}' & 0 & 0 & 0 & b_2' \\
\\
0 & 0 & a_{3,3}' & 0 & 0 & b_3' \\
\\
0 & 0 & 0 & a_{4,4}' & 0 & b_4' \\
\\
0 & 0 & 0 & 0 & a_{5,5}' & b_5' \\
\end{bmatrix}
\]
剩余系数以此化为 \(1\) 就能得到解。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1000 + 7;
const int inf = 0x3f3f3f3f;
const double eps = 1e-4;
int n;
double a[maxn][maxn];
int guass_jordan() {
for (int i = 1; i <= n; ++i) {
int maxi = i;
for (int j = i + 1; j <= n; ++j)
if (abs(a[maxi][i]) < abs(a[j][i]))
maxi = j;
if (maxi != i) swap(a[i], a[maxi]);
if (abs(a[i][i]) <= eps) return 0;
for (int j = 1; j <= n; ++j) {
if (j == i) continue;
double mul = a[j][i] / a[i][i];
for (int k = i; k <= n + 1; ++k)
a[j][k] -= mul * a[i][k];
}
}
for (int i = 1; i <= n; ++i)
a[i][n + 1] /= a[i][i];
return 1;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n + 1; ++j)
scanf("%lf", &a[i][j]);
if (guass_jordan()) {
for (int i = 1; i <= n; ++i)
printf("%.2lf\n", a[i][n + 1]);
}
else printf("No Solution\n");
return 0;
}
关于为什么系数为 \(0\) 就是无解或无限解
很显然我们可以在消元过程中遇到没有一个系数主元 \(x_i\) 就把它跳过，因为它不会影响剩下若干元的求解。
那么在高斯消元中，最后带入求解时，就会有一行的系数全是 \(0\)。此时再看等号右边，如果为 \(0\)，那就有无数解；如果不为 \(0\)，那就无解。
可以类比一元一次方程 \(kx=b\)，
\(k=0\) 时，若 \(b = 0\) 就有无数解，若 \(b \neq 0\) 就无解。
题目2
这题让具体判断是有无限解还是无解的情况。
根据上面的分析，当一个元的系数为 \(0\) 时，等号右边等于 \(0\) 就是有无数解，不等于 \(0\) 就是无解。
大体思路还是一样的，不过在寻找当前元 \(x_i\) 的系数非零方程时，不应该从第 \(i\) 个开始找，而应该是从第 \(1\) 个开始找。且如果找到的方程是在第 \(i\) 个方程前面（假如是 \(j\)），那么就要判断那个方程的主元 \(x_j\) 的系数是否为零。为零就用，不为零就不能用，因为如果不为零，主元 \(x_j\) 和 \(x_i\) 就会出现在同一个方程且系数都不为 \(0\)，此时就没有办法进行消元。
至于为什么从 \(1\) 开始，可以看看这篇题解。
注意： 无解的优先级比无数解高。即：在出现了若干行系数都为 \(0\)，而等号右边有的等于零，有有的不等于零的时，输出无解，而不是无数解。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 50 + 7;
const int inf = 0x3f3f3f3f;
const double eps = 1e-8;
int n;
double a[maxn][maxn];
int guass_jordan() {
for (int i = 1; i <= n; ++i) {
int maxi = i;
for (int j = 1; j <= n; ++j) {
if (abs(a[j][j]) > eps && j < i) continue;
if (abs(a[maxi][i]) < abs(a[j][i])) maxi = j;
}
if (maxi != i) swap(a[i], a[maxi]);
if (abs(a[i][i]) <= eps) continue;
for (int j = 1; j <= n; ++j) {
if (j == i) continue;
double mul = a[j][i] / a[i][i];
for (int k = i; k <= n + 1; ++k)
a[j][k] -= mul * a[i][k];
}
}
int res = 1;
for (int i = 1; i <= n; ++i) {
if (abs(a[i][i]) <= eps) {
if (abs(a[i][n + 1]) > eps) res = -1;
else if (res != -1) res = 0;
}
}
return res;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n + 1; ++j)
scanf("%lf", &a[i][j]);
int ans = guass_jordan();
if (ans != 1) printf("%d\n", ans);
else {
for (int i = 1; i <= n; ++i)
printf("x%d=%.2lf\n", i, a[i][n + 1] / a[i][i]);
}
return 0;
}
题目3
高斯消元求异或方程组
这种看似和解方程组无关的 “带牵连的” 数列或矩阵操作实际上也可以用高斯消元来做。
在这种问题中，我们一般会把这种 “牵连” 看作一种系数，而把数列或矩阵 【所需要进行的变化】 看成方程组等号右边的常数。
分析
就拿题目所给的第二个数据为例：
我们设 \(a[i][j]\) 表示当开关 \(j\) 受到按动时，灯泡 \(i\)（也可以说是开关 \(i\)）是否会受到牵连（\(1\) 表示会牵连，\(0\) 表示不会，至于为什么会反着记待会说）。
\[a=
\begin{bmatrix}
1 & 1 & 0 \\
\\
1 & 1 & 0 \\
\\
0 & 0 & 1 \\
\end{bmatrix}
\]
这里之所以还有一个“左上—右下”的对角线，是因为按一个开关，这个开关本身相连的灯状态也会改变。
现在让它乘以我们所要求的未知数矩阵 \(X\)（矩阵里的元素 \(x_i\) 实际上就表示是否按这个开关 \(i\)）（在这个具体问题中，两矩阵元素间的运算实际上是异或，而非加法，因题而异）。
\[a*X=
\begin{bmatrix}
1 & 1 & 0 \\
\\
1 & 1 & 0 \\
\\
0 & 0 & 1 \\
\end{bmatrix}
*
\begin{bmatrix}
x_1 \\
\\
x_2 \\
\\
x_3 \\
\end{bmatrix} =
\begin{bmatrix}
1*x_1 \bigoplus 1*x_2 \bigoplus 0*x_3 \\
\\
1*x_1 \bigoplus 1*x_2 \bigoplus 0*x_3 \\
\\
0*x_1 \bigoplus 0*x_2 \bigoplus 1*x_3 \\
\end{bmatrix}
\]
它就等于每个灯泡所要进行的变化，即：这个灯的状态是否会发生改变（这个改变包括了 由亮到暗 和 由暗到亮），将状态变化矩阵加在右边（\(1\) 表示发生改变，\(0\) 表示没发生改变）：
\[a*X=
\begin{bmatrix}
1 & 1 & 0 \\
\\
1 & 1 & 0 \\
\\
0 & 0 & 1 \\
\end{bmatrix}
*
\begin{bmatrix}
x_1 \\
\\
x_2 \\
\\
x_3 \\
\end{bmatrix} =
\begin{bmatrix}
1*x_1 \bigoplus 1*x_2 \bigoplus 0*x_3 \\
\\
1*x_1 \bigoplus 1*x_2 \bigoplus 0*x_3 \\
\\
0*x_1 \bigoplus 0*x_2 \bigoplus 1*x_3 \\
\end{bmatrix} =
\begin{bmatrix}
1 \\
\\
0 \\
\\
1 \\
\end{bmatrix}
\]
现在可以看图解释为什么要反着记录了：
总得来说，是为了方便用高斯消元的方式来求解。
很明显那一个系数矩阵，
它的纵列 \(j\) 表示：按下开关 \(j\) 哪些灯泡会受到影响；
它的横列 \(i\) 表示：第 \(i\) 个灯泡会受到哪些开关的影响。
题目所给的 “操作开关 \(i\) ”就相当于在纵列 \(i\) 通过未知数 \(x_i\) 施加一个影响，“开关 \(j\) 状态也会改变” 相当于横行 \(j\) 所代表的灯泡也会收到影响。这样符合高斯消元的格式，方便求解。
到这里分析就差不多了，只是代码实现略有不同（不过是消元时把减号运算改成了异或运算）。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 30 + 7;
const int inf = 0x3f3f3f3f;
const double eps = 1e-4;
int T;
int n;
bool a[maxn][maxn];
/*
这里的写法与上面不同：
之前用的是 "行消元"，这里用的是 "列消元"。
列消元与行消元的不同是：
行消元在遇到主元的系数为 0 时，它会直接进入下一行；
而列消元在遇到主元系数为 0 时，它依旧会停留在这一行，但进入下一列。
如此一来：
若方程组有唯一解，那么最后行数就会停留在第 n + 1 行
(为什么不是第 n 行，读者可以阅读下面代码自行理解)；
若无解或有无限解，那么最后就会停留在小于等于第 n 行的地方，
而剩下的那几行，前面的系数已经全部变成了 0，
只需要看看方程右边是否为 0，就可以判断是无解还是有无限解。
若是有无限解，剩下行的个数就是自由元的个数（自由元就是可以取任意值的未知数）。
*/
int guass() {
int r, c; // 分别代表枚举的当前行、列
for (r = 1, c = 1; c <= n; ++c) {
int maxr = r;
for (int i = r + 1; i <= n; ++i)
if (a[maxr][c] < a[i][c]) maxr = i;
if (maxr != r) swap(a[maxr], a[r]);
if (a[r][c] == 0) continue; // 当前主元系数为 0，进入下一列
for (int i = r + 1; i <= n; ++i)
for (int j = n + 1; j >= c; --j)
a[i][j] ^= a[i][c] & a[r][j];
++r; // 当前主元完成了消元，进入下一行
}
if (r <= n) {
for (int i = r; i <= n; ++i)
if (a[i][n + 1]) return 0;
}
// 由于主元只有 0 和 1 （即 "开" 和 "关"）两种取值
// 所以若自由元有 k 个，方案数就是 2^k
return 1 << (n - r + 1);
}
void sol() {
memset(a, 0, sizeof(a));
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%d", &a[i][n + 1]), a[i][i] = 1;
for (int i = 1, x; i <= n; ++i) scanf("%d", &x), a[i][n + 1] ^= x;
while (1) {
int x, y; scanf("%d%d", &x, &y);
if (x == 0 && y == 0) break;
a[y][x] = 1;
}
int ans = guass();
if (ans) printf("%d\n", ans);
else printf("Oh,it's impossible~!!\n");
}
int main() {
scanf("%d", &T);
while (T--) sol();
return 0;
}
再附上一个令笔者破防的练习题：
HDU 5755 Gambler Bo
它的题解先鸽着。
浙公网安备 33010602011771号