高斯消元
高斯消元法
高斯消元法是对多元一次方程组的通式通法，可以以 \(O(n ^ 3)\) 的时间复杂度完成对 \(n\) 元一次方程组的求解，亦可以检测其解的存在性。
高斯消元的核心在于通过加减消元来将原增广矩阵转化为上三角形式，并使用带入消元求出解向量。
约定与记法
对于一个 \(n\) 元一次方程组：
\[\left\{
\begin{matrix}
M_{1, 1} x_1 & + & M_{1, 2} x_2 & + & \cdots & + & M_{1, n} x_n & = & v_1 \\
M_{2, 1} x_2 & + & M_{2, 2} x_2 & + & \cdots & + & M_{2, n} x_n & = & v_2 \\
\vdots & & \vdots & & \ddots & & \vdots & & \vdots \\
M_{n, 1} x_n & + & M_{n, 2} x_2 & + & \cdots & + & M_{n, n} x_n & = & v_n \\
\end{matrix}
\right.
\]
我们定义其系数矩阵：
\[M = \left[
\begin{matrix}
M_{1, 1} & M_{1, 2} & \cdots & M_{1, n} \\
M_{2, 1} & M_{2, 2} & \cdots & M_{2, n} \\
\vdots & \vdots & \ddots & \vdots \\
M_{n, 1} & M_{n, 2} & \cdots & M_{n, n} \\
\end{matrix}
\right],
v = \left[
\begin{matrix}
v_1 \\
v_2 \\
\vdots \\
v_n
\end{matrix}
\right]
\]
解向量：
\[x = \left[
\begin{matrix}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{matrix}
\right]
\]
则原方程组可以简写为 \(Mx = v\)。
在高斯消元中，我们把系数矩阵 \(M\) 与常数列 \(v\) 合并在一起，写作一个 \(n\) 行 \(m\) 列的增广矩阵，其中 \(m \leftarrow n + 1\)：
\[(M|v) = \left[
\begin{array}{cccc|c}
M_{1, 1} & M_{1, 2} & \cdots & M_{1, n} & v_1 \\
M_{2, 1} & M_{2, 2} & \cdots & M_{2, n} & v_2 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
M_{n, 1} & M_{n, 2} & \cdots & M_{n, n} & v_n \\
\end{array}
\right]
\]
求解方程组
高斯消元的目标就是通过加减消元，把增广矩阵转化如下形式：
\[(M|v)' = \left[
\begin{array}{cccc|c}
1 & M_{1, 2}' & \cdots & M_{1, n}' & v_1' \\
0 & 1 & \cdots & M_{2, n}' & v_2' \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 1 & v_n' \\
\end{array}
\right]
\]
则可以自下而上地代入消元并求出解向量 \(x\)。
算法的具体流程如下：
- 取其中的一行作为主行，一般我们取系数最大的来减小误差，相关的误差分析见概率与期望。
int p = i;
for (int j = i + 1; j <= n; j++)
if (abs(M[p][i]) < abs(M[j][i])) p = j;
if (abs(M[p][i]) < eps)
printf("No Solution\n"), exit(0);
if (i != p) std::swap(M[i], M[p]);
注意 C++ 浮点数比较需要考虑舍入误差
如果该列的所有元素均为 \(0\) 则代表原方程组无解/无穷解，相关的判定会在后文介绍。
- 将主行化简为首项系数为 \(1\) 的形式。
double d = M[i][i];
for (int j = i; j <= m; j++) M[i][j] /= d;
- 用主行对矩阵剩下几行进行消元。
for (int j = i + 1; j <= n; j++)
{
d = M[j][i];
for (int k = i; k <= m; k++)
M[j][k] -= d * M[i][k];
}
- 对矩阵的 \(n\) 行进行 \(n\) 次上述操作，则矩阵可以被化为上三角形式。接下来开始带入消元求出解向量。
- 从下而上逆序求出解向量。
for (int i = n; i; i--)
{
x[i] = M[i][m];
for (int j = i + 1; j <= n; j++)
x[i] -= M[i][j] * x[j];
}
则原方程组求解完毕。
模板题：
#include <bits/extc++.h>
#define inline __always_inline
template <typename T> inline void read(T &x)
{
char ch; bool flag = false;
for (ch = getchar(); !isdigit(ch); ch = getchar()) flag = ch == '-';
for (x = 0; isdigit(ch); ch = getchar()) x = x * 10 + ch - '0';
if (flag) x = -x;
}
const int MaxN = 105;
const double eps = 1e-9;
int n, m; double M[MaxN][MaxN], x[MaxN];
int main()
{
read(n), m = n + 1;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= m; j++)
read(M[i][j]);
for (int i = 1; i <= n; i++)
{
int p = i;
for (int j = i + 1; j <= n; j++)
if (abs(M[p][i]) < abs(M[j][i])) p = j;
if (abs(M[p][i]) < eps)
printf("No Solution\n"), exit(0);
if (i != p) std::swap(M[i], M[p]);
double d = M[i][i];
for (int j = i; j <= m; j++) M[i][j] /= d;
for (int j = i + 1; j <= n; j++)
{
d = M[j][i];
for (int k = i; k <= m; k++)
M[j][k] -= d * M[i][k];
}
}
for (int i = n; i; i--)
{
x[i] = M[i][m];
for (int j = i + 1; j <= n; j++)
x[i] -= M[i][j] * x[j];
}
for (int i = 1; i <= n; i++) printf("%.2lf\n", x[i]);
return 0;
}
解的存在性判定
在上述的消元过程中，如果在一次消元中，当前列的所有系数均为 \(0\)，则说明方程组可能无解/无穷解。
可以对原算法做如下修改来实现对解的存在性的具体判定：
使用一个指针 top 来指向当前修改的行。
int top = 1;
---
int p = top;
for (int j = top + 1; j <= n; j++)
if (abs(M[p][i]) < abs(M[j][i])) p = j;
if (abs(M[p][i]) < eps) continue; // 无解/无穷解，跳过这一列继续消元。
if (top != p) std::swap(M[top], M[p]);
double d = M[top][i];
for (int j = i; j <= m; j++) M[top][j] /= d;
for (int j = top + 1; j <= n; j++)
{
d = M[j][i];
for (int k = i; k <= m; k++)
M[j][k] -= d * M[top][k];
}
top++;
若当前这一列系数均为 \(0\)，则继续消元下一列，直至完全消元至上三角形式。
- 若无解，则消元完成后 \(\exists i, M_{i, i} = 0 \wedge M_{i, m} \neq 0\)
- 若无穷解，则消元完成后 \(\forall i, M_{i, i} = 0 \wedge M_{i, m} = 0\)
#include <bits/extc++.h>
#define inline __always_inline
template <typename T> inline void read(T &x)
{
char ch; bool flag = false;
for (ch = getchar(); !isdigit(ch); ch = getchar()) flag = ch == '-';
for (x = 0; isdigit(ch); ch = getchar()) x = x * 10 + ch - '0';
if (flag) x = -x;
}
const int MaxN = 55;
const double eps = 1e-9;
int n, m; double M[MaxN][MaxN], v[MaxN];
void debug()
{
for (int i = 1; i <= n; i++, putchar('\n'))
for (int j = 1; j <= m; j++)
printf("%.2lf ", M[i][j]);
}
int main()
{
read(n), m = n + 1;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= m; j++)
read(M[i][j]);
int top = 1;
for (int i = 1; i <= n; i++)
{
int p = top;
for (int j = top + 1; j <= n; j++)
if (abs(M[p][i]) < abs(M[j][i])) p = j;
if (abs(M[p][i]) < eps) continue;
if (top != p) std::swap(M[top], M[p]);
double d = M[top][i];
for (int j = i; j <= m; j++) M[top][j] /= d;
for (int j = top + 1; j <= n; j++)
{
d = M[j][i];
for (int k = i; k <= m; k++)
M[j][k] -= d * M[top][k];
}
top++;
}
if (top <= n)
{
for (int i = top; i <= n; i++)
if (abs(M[i][m]) > eps) printf("-1"), exit(0);
printf("0"), exit(0);
}
for (int i = n; i; i--)
{
v[i] = M[i][m];
for (int j = i + 1; j <= n; j++)
v[i] -= M[i][j] * v[j];
}
for (int i = 1; i <= n; i++) printf("x%d=%.2lf\n", i, v[i]);
return 0;
}
浙公网安备 33010602011771号