用高斯消元法求解多元线性方程组
多元线性方程组
形式为：
\[\begin{cases}
a_{11}x_1+a_{12}x_2+\dots+a_{1n}x_n=b_1 \\
a_{21}x_1+a_{22}x_2+\dots+a_{2n}x_n=b_2 \\
\vdots \\
a_{n1}x_1+a_{n2}x_2+\dots+a_{nn}x_n=b_n \\
\end{cases}
\]
求解 \(x_1,x_2,\dots,x_n\).
解的可能一共有三种：
\[\begin{cases}
1.有唯一解 \\
2.无解 \\
3.有无穷多组解 \\
\end{cases}
\]
解法：高斯消元法
初等行列变换操作：
（以下操作是对原方程组的等价变形）
把某一行乘一个非零的数；（等式两边同时乘一个非0的数）
交换某两行；（交换某两个方程的位置）
把某行的若干倍加到另一行上去。（某个方程乘若干倍加到另一个方程上去）
步骤：
枚举每一列：
找到剩下行中当前这一列绝对值最大的一行；
操作2：将这一行换到最上面去；
操作1：将该行的第一个数变成1，即每个系数都等比缩小（或放大）；
操作3：用这一行将下面所有行的当前列消成0。
举例：
1.有唯一解
\[\begin{cases}
-x_1-x_2+2x_3=7 \\
2x_1+x_2-3x_3=-9 \\
x_1+2x_2-x_3=-6 \\
\end{cases}
\]
转化为矩阵：
-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\) 会枚举完每一行，最终矩阵会变成上面的三角形的形式。也就是将方程组化简为这样：
\[\begin{cases}
x_1+0.5x_2-1.5x_3=-4.5 \\
\qquad \qquad x_2+\frac{1}{3}x_3=-1 \\
\qquad \qquad \qquad \quad x_3=3 \\
\end{cases}
\]
我们就可以倒着一步步推出每个解。最终每个解就是矩阵每行的最后一个数字，即 \(x_i=a[i][n+1]\)。
2.无解
\[\begin{cases}
x_1=0 \\
x_3=1 \\
x_3=0 \\
\end{cases}
\]
其对应的矩阵表示为：
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
也就等价于：
\[\begin{cases}
x_1=0 \\
x_3=1 \\
0=-1 \\
\end{cases}
\]
由此我们看出，经过等价变形，我们得出了一个显然不成立的式子 \(0=-1\)，所以这个方程组无解。
也就是说，当出现 \(0=非零数\) 这样形式的式子的时候，就是方程组无解的时候。
3.有无穷组解
\[\begin{cases}
x_1+x_3=2 \\
x_1+x_2=3 \\
2x_1+2x_2=6 \\
\end{cases}
\]
首先我们会发现，(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
#include
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.解异或线性方程组
形式为：
\[\begin{cases}
a_{11}x_1 \land a_{12}x_2\land \dots \land a_{1n}x_n=b_1 \\
a_{21}x_1 \land a_{22}x_2 \land \dots \land a_{2n}x_n=b_2 \\
\vdots \\
a_{n1}x_1 \land a_{n2}x_2 \land \dots \land a_{nn}x_n=b_n \\
\end{cases}
\]
其中 ^ 表示异或运算（XOR），方程组的系数 \(a\) 和常数 \(b\) 均为 \(0\) 或 \(1\)，未知数可能的值也为 \(0\) 或 \(1\)。
求解： \(x_1,x_2,\dots,x_n\).
思路
与一般的高斯消元法相同，只是稍作改动。
不再需要 double类型，使用 bool 类型即可；
不再需要等比例缩小至系数为 \(1\)，因为系数不是 \(0\) 就是 \(1\)；
加减消元可以用 ^ 运算来代替，因为异或运算俗称不进位的加法运算。
乘法运算可以用 & 运算来代替，众所周知位运算是最快的运算。
举例
\[\begin{cases}
x_1 \land x_2=1 \\
x_2 \land x_3=0 \\
x_1=1 \\
\end{cases}
\]
表示为：
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
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;
}