[算法] 高斯消元及其应用

主要问题

对于一个线性的方程组求解。

假设这个方程有 \(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 <cmath>
#include <cstdio>
#include <algorithm>
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 <cmath>
#include <cstdio>
#include <algorithm>
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;
}