1 基础知识
高斯消元法,用来求解线性方程组的解,
{
a
11
x
1
+
a
12
x
2
+
⋯
+
a
1
n
x
n
=
b
1
a
21
x
1
+
a
22
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
\left \{ \begin{matrix} a_{11}x_1 + a_{12}x_2+\cdots+a_{1n}x_n=b_1 \\ a_{21}x_1 + a_{22}x_2+\cdots+a_{2n}x_n=b_2 \\ \cdots \\ a_{n1}x_1 + a_{n2}x_2+\cdots+a_{nn}x_n=b_n \\ \end{matrix} \right.
⎩
⎨
⎧a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋯an1x1+an2x2+⋯+annxn=bn
初等行列变换:
- 某行乘以k倍,即row[i] = k * row[i]
- 交换某两行,即row[i], row[j] = row[j], row[i]
- 某行加上另一行的k倍,即row[i] = row[i] + k * row[j]
高斯消元法代码实现的关键步骤:
- 枚举每一列c,进行以下操作:
– 找到绝对值最大的一行(注意,要从不固定的行中寻找)。
– 将该行换到最上面。
– 将该行的第1个数变成1。
– 将下面所有行的第c列消成0。 - 经过上述步骤,行列式已转变为下三角行列式,需要从第n行(或者第n-1行)开始把未知数 x i x_i xi的值代入方程进行求解,依次求解出 x n x_n xn、 x n − 1 x_{n-1} xn−1、……、 x 1 x_1 x1的值。
将上述操作转换为代码,如下,
#include <iostream>
#include <cmath>
using namespace std;
const int N = 110;
const double eps = 1e-6;
double a[N][N];
int n;
int gauss() {
//从0到n-1遍历每一列,把行列式消成下三角形式
int r, c;
for (r = 0, c = 0; c < n; ++c) {
//正在处理第r行第c列,a[r][c]
//step1:找到r~n-1行中第c列绝对值最大的行,记为t
int t = r;
for (int i = r + 1; i < n; ++i) {
if (fabs(a[i][c]) > fabs(a[t][c])) {
t = i;
}
}
//step2:交换第r行和第t行的第c~n列的元素值
for (int j = c; j <= n; ++j) {
swap(a[r][j], a[t][j]);
}
//step3:将第r行的变量x_c前面的系数变为1,即将a[r][c]变为1
if (fabs(a[r][c]) < eps) {
continue;
}
for (int j = n; j >= c; --j) {
a[r][j] = a[r][j] / a[r][c];
}
//step4:将第r+1~n-1行中的第c列给去掉,即将第r行的某倍数加到下面的行中
for (int i = r + 1; i < n; ++i) {
//将第i行中的第c列给去掉(即x_c的系数为0)
for (int j = n; j >= c; --j) {
a[i][j] -= a[i][c] * a[r][j]; //将第r行的-a[i][c]倍加到第i行
}
}
r += 1; //有效方程数加1
}
if (r < n) {
for (int i = r; i < n; ++i) {
if (fabs(a[i][n]) > eps) {
return 2; //0等于非0,无解。
}
}
return 1; //有效方程数小于未知量的个数,无穷多组解。
}
//step5:逆向求解x_{n-1},x_{n-2}...x_0
for (int i = n - 1; i >= 0; --i) {
//求解未知数x_i,第i行第i列
for (int j = i + 1; j < n; ++j) {
a[i][n] -= a[i][j] * a[j][n]; //a[j][n]表示x_j的值
}
}
return 0;
}
int main() {
cin >> n;
//总共n行、n+1列
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n + 1; ++j) {
cin >> a[i][j];
}
}
int ans = gauss();
if (ans == 0) {
for (int i = 0; i < n; ++i) {
printf("%.2lf\n", a[i][n]); //xi的解为a[i][n]
}
} else if (ans == 1) {
cout << "Infinite group solutions" << endl;
} else {
cout << "No solution" << endl;
}
return 0;
}
2 模板
// a[N][N]是增广矩阵
int gauss()
{
int c, 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) continue;
for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]); // 将绝对值最大的行换到最顶端
for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c]; // 将当前行的首位变成1
for (int i = r + 1; i < n; i ++ ) // 用当前行将下面所有的列消成0
if (fabs(a[i][c]) > eps)
for (int j = n; j >= c; j -- )
a[i][j] -= a[r][j] * a[i][c];
r ++ ;
}
if (r < n)
{
for (int i = r; i < n; i ++ )
if (fabs(a[i][n]) > eps)
return 2; // 无解
return 1; // 有无穷多组解
}
for (int i = n - 1; i >= 0; i -- )
for (int j = i + 1; j < n; j ++ )
a[i][n] -= a[i][j] * a[j][n];
return 0; // 有唯一解
}
3 工程化
题目1:高斯消元法求解异或线性方程组。
#include <iostream>
using namespace std;
const int N = 110;
int n;
int a[N][N];
int gauss() {
int r, c;
for (r = c = 0; c < n; ++c) {//枚举列
//找到r~n行中a[][c]不为零的行
int t = r;
for (int i = r; i < n; ++i) {
if (a[i][c]) {
t = i;
break;
}
}
if (!a[t][c]) continue;
//交换第r行和第t行,注意j>=c
for (int j = c; j <= n; ++j) swap(a[r][j], a[t][j]);
//将第r行之下的第c列均消成零
for (int i = r + 1; i < n; ++i) {
//将第i行的第c列均消成零
if (a[i][c]) {
//将第r行异或下来即可
for (int j = c; j <= n; ++j) {
a[i][j] ^= a[r][j];
}
}
}
r++;
}
if (r < n) {
for (int i = r; i < n; ++i) {
if (a[i][n]) {
return 2;
}
}
return 1;
}
//从第n-1行开始把答案计算出来,存储在a[i][n]
for (int i = n - 1; i >= 0; --i) {
//计算x[i],存储在a[i][n]中
for (int j = i + 1; j < n; ++j) {
a[i][n] = a[i][n] ^ a[i][j] & a[j][n];
}
}
return 0;
}
int main() {
cin >> n;
for (int i = 0; i < n; ++i) {
for (int j = 0; j <= n; ++j) {
cin >> a[i][j];
}
}
int res = gauss();
if (res == 0) {
for (int i = 0; i < n; ++i) {
cout << a[i][n] << endl;
}
} else if (res == 1) {
cout << "Multiple sets of solutions" << endl;
} else {
cout << "No solution" << endl;
}
return 0;
}