高斯消元法求解多元线性方程组
用高斯消元法求解多元线性方程组
多元线性方程组
形式为：
求解 \(x_1,x_2,\dots,x_n\).
解的可能一共有三种：
解法：高斯消元法
初等行列变换操作：
（以下操作是对原方程组的等价变形）
- 把某一行乘一个非零的数；（等式两边同时乘一个非0的数）
- 交换某两行；（交换某两个方程的位置）
- 把某行的若干倍加到另一行上去。（某个方程乘若干倍加到另一个方程上去）
步骤：
枚举每一列：
- 找到剩下行中当前这一列绝对值最大的一行；
- 操作2：将这一行换到最上面去；
- 操作1：将该行的第一个数变成1，即每个系数都等比缩小（或放大）；
- 操作3：用这一行将下面所有行的当前列消成0。
举例：
1.有唯一解
转化为矩阵：
| -1 | -1 | 2 | 7 |
|---|---|---|---|
| 2 | 1 | -3 | -9 |
| 1 | 2 | -1 | -6 |
对于第 1 列：
- 先找到 \(2\) 最大；
- 将第 \(2\) 行换到第 \(1\) 行；
| 2 | 1 | -3 | -9 |
|---|---|---|---|
| -1 | -1 | 2 | 7 |
| 1 | 2 | -1 | -6 |
- 使第一个数变成 \(1\)，同时除以 \(2\)；
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| -1 | -1 | 2 | 7 |
| 1 | 2 | -1 | -6 |
- 用这一行将下面所有行的当前列消成 \(0\)：第 \(2\) 行减第 \(1\) 行的 \(-1\) 倍；第3行减第1行的 \(1\) 倍。
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| 0 | -0.5 | 0.5 | 2.5 |
| 0 | 1.5 | 0.5 | -1.5 |
对于第 2 列：
- 从第 2 行开始找，第 2 列绝对值最大的是1.5；
- 将第 \(3\) 行换到第 2 行；
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| 0 | 1.5 | 0.5 | -1.5 |
| 0 | -0.5 | 0.5 | 2.5 |
- 使第一个数变成 \(1\)，同时除以 \(1.5\)；
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| 0 | 1 | 1/3 | -1 |
| 0 | -0.5 | 0.5 | 2.5 |
- 用这一行将下面所有行的当前列消成 \(0\)：第 \(3\) 行减第 \(1\) 行的 \(-0.5\) 倍。
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| 0 | 1 | 1/3 | -1 |
| 0 | 0 | 2/3 | 2 |
对于第 3 列：
- 从第 3 行开始找，第 3 列绝对值最大 的是 \(\frac{2}{3}\)；
- 将第 3 行换到第 3 行，相当于不变；
- 使第 3 行第一个数变成 \(1\)，同时除以 \(\frac{2}{3}\)；
| 1 | 0.5 | -1.5 | -4.5 |
|---|---|---|---|
| 0 | 1 | 1/3 | -1 |
| 0 | 0 | 1 | 3 |
- 没有下一行，跳过此步。
所有有唯一解的方程组，经过此 \(4\) 步的循环操作之后，行数 \(r\) 会枚举完每一行，最终矩阵会变成上面的三角形的形式。也就是将方程组化简为这样：
我们就可以倒着一步步推出每个解。最终每个解就是矩阵每行的最后一个数字，即 \(x_i=a[i][n+1]\)。
2.无解
其对应的矩阵表示为：
| 1 | 0 | 0 | 0 |
|---|---|---|---|
| 0 | 0 | 1 | 1 |
| 0 | 0 | 1 | 0 |
枚举到第 2 列时，从第 2 行看起，会发现绝对值最大的是0；
于是直接从第 3 列的第 2 行看起，找到绝对值最大的第 2 行的 1，换到第 2 行不必交换了。然后把下面的全部变成 \(0\)，也就是用第 3 行减去第 2 行。结果为：
| 1 | 0 | 0 | 0 |
|---|---|---|---|
| 0 | 0 | 1 | 1 |
| 0 | 0 | 0 | -1 |
也就等价于：
由此我们看出，经过等价变形，我们得出了一个显然不成立的式子 \(0=-1\)，所以这个方程组无解。
也就是说，当出现 \(0=非零数\) 这样形式的式子的时候，就是方程组无解的时候。
3.有无穷组解
首先我们会发现，(3)式是由(2)式乘2得到，相当于这个方程组只有两个式子，却有3个未知数，显然有无穷组解。
我们看看程序运行过程是怎样的：
| 1 | 0 | 1 | 2 |
|---|---|---|---|
| 1 | 1 | 0 | 3 |
| 2 | 2 | 0 | 6 |
对于第 1 列：
首先交换1、3行，将 \(x_1\) 的系数变为 \(1\)，变成： \(\begin{cases} x_1+x_2=3 \\ x_1+x_2=3 \\ x_1+x_3=2 \\ \end{cases}\)
然后将2、3行的 \(x_1\) 项消掉，即(2)式减(1)式，(3)式减((1)式，结果为：\(\begin{cases} x_1+x_2=3 \\ 0=0 \\ -x_2+x_3=-1 \\ \end{cases}\)，即：
| 1 | 1 | 0 | 3 |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 0 | -1 | 1 | -1 |
对于第 2 列：
首先交换2、3行，将 \(x_2\) 的系数变为 \(1\)，变成：\(\begin{cases} x_1+x_2=3 \\ x_2-x_3=1 \\ 0=0 \\ \end{cases}\)；下面的 \(x_2\) 项系数全部为 \(0\)，因此无需改动。即：
| 1 | 1 | 0 | 3 |
|---|---|---|---|
| 0 | 1 | -1 | 1 |
| 0 | 0 | 0 | 0 |
枚举的行数 \(r\) 没有枚举到最后一行就已经完毕，因此属于非正常情况，又没有出现 0=非零数 这样的无解情况，因此是无穷组解。
代码
#include <iostream>
#include <cmath>
using namespace std;
const int N = 105;
const double eps = 1e-6;
int n;
double a[N][N];
int gauss() { //高斯消元法
int r, c, i, j; //r为行，c为列
for (r = 1, c = 1; c <= n; c++) {
int t = r;
for (i = r; i <= n; i++) //循环找到当前列绝对值最大的行t
if (fabs(a[i][c]) > fabs(a[t][c])){
t = i;
break;
}
if (fabs(a[t][c]) < eps) //为0
continue;
if (t != r) //有更大的才交换
for (i = c; i <= n + 1; i++) //交换到当前行，从c列开始交换是因为前面都是0
swap(a[r][i], a[t][i]);
if (a[r][c] != 1) //不是1才需要变成1
for (i = n + 1; i >= c; i--) //将第c列的数字变为1，第一个数字要最后变
a[r][i] /= a[r][c];
for (i = r + 1; i <= n; i++) //将下面每行的第c列变成0
if (fabs(a[i][c]) > eps) //不是0的才变，是0就不变了
for (j = n + 1; j >= c; j--)
a[i][j] -= a[r][j] * a[i][c];
r++;
}
if (r <= n) { //说明不是正常情况，正常r=n+1
for (i = r; i <= n; i++)
if (fabs(a[i][n + 1]) > eps)
return -1; //无解
return 2; //无穷组解
}
for (i = n; i > 0; i--) //倒着一步步推出每个解
for (j = i + 1; j <= n; j++)
a[i][n + 1] -= a[i][j] * a[j][n + 1];
return 1; //有唯一解
}
int main() {
int i, j;
scanf("%d", &n);
for (i = 1; i <= n; i++)
for (j = 1; j <= n + 1; j++) //输入的是一个n*(n+1)的矩阵
scanf("%lf", &a[i][j]);
int t = gauss();
if (t == -1) { //无解
puts("No solution");
} else if (t == 2) { //有无穷组解
puts("Infinite group solutions");
} else { //有唯一解
for (i = 1; i <= n; i++)
printf("%.2lf\n", a[i][n + 1]);
}
return 0;
}
变式
1.解异或线性方程组
形式为：
其中 ^ 表示异或运算（XOR），方程组的系数 \(a\) 和常数 \(b\) 均为 \(0\) 或 \(1\)，未知数可能的值也为 \(0\) 或 \(1\)。
求解： \(x_1,x_2,\dots,x_n\).
思路
与一般的高斯消元法相同，只是稍作改动。
- 不再需要
double类型，使用bool类型即可； - 不再需要等比例缩小至系数为 \(1\)，因为系数不是 \(0\) 就是 \(1\)；
- 加减消元可以用
^运算来代替，因为异或运算俗称不进位的加法运算。 - 乘法运算可以用
&运算来代替，众所周知位运算是最快的运算。
举例
表示为：
| 1 | 1 | 0 | 1 |
|---|---|---|---|
| 0 | 1 | 1 | 0 |
| 1 | 0 | 0 | 1 |
将第 \(3\) 行的第 \(1\) 列变成 \(0\)，就是将第 \(1\) 行与第 \(3\) 行的系数做异或运算再放入第 \(3\) 行，变成：\(\begin{cases} x_1 \land x_2=1 \\ x_2 \land x_3=0 \\ x_2=0 \\ \end{cases}\)
| 1 | 1 | 0 | 1 |
|---|---|---|---|
| 0 | 1 | 1 | 0 |
| 0 | 1 | 0 | 0 |
将第 \(3\) 行的第 \(2\) 列变成 \(0\)，就是将第 \(2\) 行与第 \(3\) 行的系数做异或运算再放入第 \(3\) 行，变成：\(\begin{cases} x_1 \land x_2=1 \\ x_2 \land x_3=0 \\ x_3=0 \\ \end{cases}\)
| 1 | 1 | 0 | 1 |
|---|---|---|---|
| 0 | 1 | 1 | 0 |
| 0 | 0 | 1 | 0 |
最后倒推回去即可。这里注意：异或运算的逆运算还是异或运算。因此要求 \(x_2\)，只需求 \(x_2=0 \land x_3=0 \land 0=0\) 即可。
最终可得到，此方程组的解为：\(\begin{cases} x_1=1 \\ x_2=0 \\ x_3=0 \\ \end{cases}\)
代码
#include <iostream>
using namespace std;
const int N = 105;
int n;
bool a[N][N];
int gauss() { //高斯消元法
int r, c, i, j; //r为行，c为列
for (r = 1, c = 1; c <= n; c++) {
int t = r;
for (i = r; i <= n; i++) //循环找到当前列绝对值最大的行t
if (a[i][c]){ //是1就可以
t = i;
break;
}
if (!a[t][c]) continue; //为0
if (t != r) //有更大的才交换
for (i = c; i <= n + 1; i++) //交换到当前行，从c列开始交换是因为前面都是0
swap(a[r][i], a[t][i]);
for (i = r + 1; i <= n; i++) //将下面每行的第c列变成0
if (a[i][c]) //是1才需要变
for (j = n + 1; j >= c; j--)
a[i][j] ^= a[r][j];
r++;
}
if (r <= n) { //说明不是正常情况，正常r=n+1
for (i = r; i <= n; i++)
if (a[i][n + 1])
return -1; //无解
return 2; //多组解
}
for (i = n; i > 0; i--) //倒着一步步推出每个解
for (j = i + 1; j <= n; j++)
a[i][n + 1] ^= a[i][j] & a[j][n + 1];
return 1; //有唯一解
}
int main() {
int i, j;
scanf("%d", &n);
for (i = 1; i <= n; i++)
for (j = 1; j <= n + 1; j++) //输入的是一个n*(n+1)的矩阵
cin >> a[i][j]; //此处注意，如果想用scanf读入，就不能使用bool数组了，应换成int
int t = gauss();
if (t == -1) { //无解
puts("No solution");
} else if (t == 2) { //有多组解
puts("Multiple sets of solutions");
} else { //有唯一解
for (i = 1; i <= n; i++)
printf("%d\n", a[i][n + 1]);
}
return 0;
}
浙公网安备 33010602011771号