[算法] 高斯消元及其应用
主要问题
对于一个线性的方程组求解。
假设这个方程有 \(n\) 个，则时间复杂度为 \(O(n^3)\)。
有些题目的 \(dp\) 状态有后效性，但是对于线性的方程，可以用高斯消元进行计算。
解决方法
高斯消元法的思路是：通过消元运算，直到得到一个只关于 \(x_1\) 的式子，只关于 \(x_1,x_2\) 的式子，只关于 \(x_1,x_2,x_3\) 的式子，然后回代计算即可。
如下：最初有一个 \(n\) 元线性方程
\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1 \\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n=b_2 \\ \vdots \\ a_{n,1}x_1+a_{n,2}x_2+...+a_{n,n}x_n=b_n \\ \end{cases} \]将上述方程转换为
\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1 \\ a_{1,1}x_1+\frac{a_{2,2}a_{1,1}}{a_{2,1}}x_2+...+\frac{a_{2,n}a_{1,1}}{a_{2,1}}x_n=b_2 \\ \vdots \\ a_{1,1}x_1+\frac{a_{n,2}a_{1,1}}{a_{n,1}}x_2+...+\frac{a_{n,n}a_{1,1}}{a_{n,1}}x_n=b_n \\ \end{cases} \]最后通过减法消去后 \(n-1\) 个方程的 \(x_1\)：
\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1 \\ (\frac{a_{2,2}a_{1,1}}{a_{2,1}}-a_{1,1})x_2+...+(\frac{a_{2,n}a_{1,1}}{a_{2,1}}-a_{1,n})x_n=b_2-b_1 \\ \vdots \\ (\frac{a_{n,2}a_{1,1}}{a_{n,1}}-a_{1,1})x_2+...+(\frac{a_{n,n}a_{1,1}}{a_{n,1}}-a_{1,n})x_n=b_n-b_1 \\ \end{cases} \]如此，对于剩下的 \(n-1\) 个方程又是一个新的线性方程求解问题，继续递归求解，最后回代。
若对于当前的子问题中，有一次计算中的 \(x_n\) 就已经消失了，那么原问题中的这一未知数没有唯一解。
若当前中存在一个方程的各项系数都为 \(0\) ，且该条方程的常数项不为 \(0\) ，则一定无解。
当然，代码具体实现的时候并不需要写递归，直接用循环模拟即可。
实现中还有一个问题，每次求解的 \(a_{1,1}\) 是否会对答案产生影响？
可以发现，若 \(a_{1,1}\) 太小，那么会导致精度问题，所以每次需要选取主元素：最大的 \(a_{i,1}\) 所有方程以它为基础来进行消元。
- Code:
#include
#include
#include
using namespace std;
#define eps 1e-7
void Read(int &n) {
n = 0; bool f = 0; char c = getchar();
while (c < '0' || c > '9') {
if (c == '-') f = 1;
c = getchar();
}
while (c >= '0' && c <= '9') {
n = (n << 1) + (n << 3) + (c ^ 48);
c = getchar();
}
if (f) n = -n;
}
const int MAXN = 1e2 + 5;
double a[MAXN][MAXN], ans[MAXN];
int n;
int main() {
Read(n);
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n + 1; j++) scanf("%lf", &a[i][j]);
for (int i = 1; i <= n; i++) {
int r = i;
for (int j = i + 1; j <= n; j++)
if (fabs(a[r][i]) < fabs(a[j][i])) r = j;//选取主元素
if (fabs(a[r][i]) < eps) return printf("No Solution"), 0;//不存在唯一解
if (i != r) swap(a[i], a[r]);
double tmp = a[i][i];//带入消元
for (int j = i; j <= n + 1; j++) a[i][j] /= tmp;
for (int j = i + 1; j <= n; j++) {
double tmp = a[j][i];
for (int k = i; k <= n + 1; k++) a[j][k] -= a[i][k] * tmp;
}
}
ans[n] = a[n][n + 1];//回代求解
for (int i = n - 1; i >= 1; i--) {
ans[i] = a[i][n + 1];
for (int j = i + 1; j <= n; j++) ans[i] -= a[i][j] * ans[j];
}
for (int i = 1; i <= n; i++) printf("%.2lf\n", ans[i]);
return 0;
}
模板题目链接。
[JSOI2008]球形空间产生器
题目大意
在 \(n\) 维空间中，给出 \(n+1\) 个点，求出圆心的坐标。
规定：
- 球心：到球面上任意一点距离都相等的点。
- 距离：设两个 \(n\) 维空间上的点 \(A, B\) 的坐标为 \((a_1, a_2, \cdots , a_n), (b_1, b_2, \cdots , b_n)\)，则AB的距离定义为：\(dist = \sqrt{ (a_1-b_1)^2 + (a_2-b_2)^2 + \cdots + (a_n-b_n)^2 }\)
思路
设坐标数组为 \(a\)，由题意得：
\(\sum_{j=0}^{n} (a_{i,j}-x_j)^2=C\)
化简得：
\(\sum_{j=1}^{n}(a_{i,j}^2-a_{i+1,j}^2-2x_j(a_{i,j}-a_{i+1,j}))=0\)
将常数项与未知数分离：
\(\sum_{j=1}^{n}2(a_{i,j}-a_{i+1,j})x_j=\sum_{j=1}^{n}(a_{i,j}^2-a_{i+1,j}^2)\)
等价转化为了 \(n\) 个线性方程，求解即可。
Code
#include
#include
#include
using namespace std;
#define eps 1e-7
const int MAXN = 1e2 + 5;
double a[MAXN][MAXN], b[MAXN][MAXN], ans[MAXN];
int n;
int main() {
scanf("%d", &n);
for (int i = 1; i <= n + 1; i++)
for (int j = 1; j <= n; j++) scanf("%lf", &b[i][j]);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
a[i][j] = 2 * (b[i][j] - b[i + 1][j]);
a[i][n + 1] += b[i][j] * b[i][j] - b[i + 1][j] * b[i + 1][j];
}
}
for (int i = 1; i <= n; i++) {
int r = i;
for (int j = i + 1; j <= n; j++)
if (fabs(a[r][i]) < fabs(a[j][i])) r = j;
if (i != r) for (int j = 1; j <= n + 1; j++) swap(a[i][j], a[r][j]);
double tmp = a[i][i];
for (int j = i; j <= n + 1; j++) a[i][j] /= tmp;
for (int j = i + 1; j <= n; j++) {
double tmp = a[j][i];
for (int k = i; k <= n + 1; k++) a[j][k] -= a[i][k] * tmp;
}
}
ans[n] = a[n][n + 1];
for (int i = n - 1; i >= 1; i--) {
ans[i] = a[i][n + 1];
for (int j = i + 1; j <= n; j++) ans[i] -= a[i][j] * ans[j];
}
for (int i = 1; i <= n; i++) printf("%.3lf ", ans[i]);
return 0;
}