例题:[NOIP2004 提高组] 虫食算
在这道题中,我们要解一个多元一次方程组,即对于
这样的数式计算,要找到所有字母对应的数值,使对于这样的 进制加法成立。每个字母对应一个数字。
很容易可以转换为方程组。问题就在于,如何解多元方程组?
答案就是,高斯消元。
做法
首先,可以把一个多元方程转换为一个矩阵:
想象平时解一个多元方程时使用的方法:消元。
那么如何对一个多元方程进行消元呢?
高斯消元选用加减消元法。
想象一个方程:
,
于是方程变为
通过这样的变形,方程组将变成一个好看的形状——斜三角形 /kk。
观察到,在这种形状中,第 个方程中,所有 的 系数都为 。
这样的形状有什么用呢?
最下方的式子形式为 ,可以获得一组解;
通过代入的方式,倒数第二行的方程 由于已经知道了 的解,可以直接带入求 的值。
通过带入和消元,最终将会把所有方程全部消干净,于是就能得到多元方程的解。
将方程化为矩阵,上面的例子可以解释为这样:
对矩阵进行初等变换,
可以变成这样子
只要将矩阵化为左下角全为 就肯定有解了。
那么,问题就在于,如何把这一并非非常客观的过程转化为代码实现呢?
实现
要得到一般规律,我们需要寻找一个普遍的实现,来对这个实现进行总结。
观察第一次变形。
由于倒三角形式(简化阶梯型方程组)的要求,必然有 行及之后的所有第 列之前都是 。
比如说当 时,第二行及之后的所有行的第一列要求变成 。
那么如何做到这一点呢?
每次按顺序(从 开始)依次选中 ,对于所有 行及之后的所有行进行加减消元。
由于是依次进行消元的,所以选择第 行进行加减消元,理由是刚好比下面所有行都多一个系数(因为 ta 在前面消元的过程中已经消掉了 及之前的所有系数。如:
现在已经完成了第一轮消掉 的过程了。
那么如何对于第三第四行依次消掉 的呢?
选中第二行,刚好比下面都多要求保留一个系数。(在这一次中第二行不变,因为已经满足条件了)
必然可以通过乘上系数的方式使得 (这样就能通过做减法消掉系数)。
就如上面的例子,第三行乘 ,就变成了这个样子:
于是做减法。
为什么会出现小数呢(反而上面的示例没有)?
因为上面是手算——为了避免分数而调换了行之间的顺序(其实差不多,对于电脑来说,只不过是浮点而已)。
况且有些时候解必然是浮点数。
第四行同理,乘上 ,
变成
做减法,
那么将会又消掉两个数!
不断的进行下去,就可以获得解!
还有一个问题。如果在消元的过程中出现上一行的第 个数字刚刚好为 ,那么很不走运的,这一个方程并不能用来做消元。
那么该如何调换各个方程之间的顺序,使得每一次操作的时候都有一个方程可以去加减消除掉?
可以在对每列做消元的之前做个预处理。
对于每一个 ,都要保证 。
如果出现 ,那么就需要找一个 。
于是可以寻找 。
(为什么这里 下界 为 呢?因为不能够动前面的方程啊,不然前面就不满足条件了)
这时候就有聪明的小脑瓜要问了:
如果找不到呢?
那就说明这个元 没有定值呗。
那就有更聪明的小脑瓜又要问了:
如果在寻找之前,已经有一个方程中有确定的 的系数被前面的操作调整到了数列的前面(也就是说在这个时候用不了了),到后面的时候发现没有系数为 的方程了,那这样又是什么情况呢?
实际上 这玩意我也没搞太明白 XDD 反正正确性显然(
那么,构造做法到此结束了,下面就要开始实现了。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
| #include <cstdio> #include <cstring> #include <cmath> #include <vector> #include <algorithm>
#define gc getchar() template<typename T> void read(T &r) { r = 0; static char ch, last; ch = gc, last = 'z'; while (ch < '0' || ch > '9') last = ch, ch = gc; while (ch >= '0' && ch <= '9') r = (r << 1) + (r << 3) + (ch ^ 48), ch = gc; r = (last == '-') ? -r : r; }
template<typename T, typename...Ts> void read(T &arg, Ts&...arg_left) { read(arg); read(arg_left...); }
double mat[202][202]; double ans[202];
int main() { int n; read(n);
for (int i = 0; i < n; i++) { for (int j = 0; j < n + 1; j++) { scanf("%lf", &mat[i][j]); } }
for (int k = 0; k < n; k++) { bool flag = false; for (int j = k; j < n; j++) { if (mat[j][k] != 0) { std::swap(mat[k], mat[j]); flag = true; break; } } if (!flag) { puts("No Solution"); return 0; } for (int i = k + 1; i < n; i++) { if (fabs(mat[i][k]) <= 1e-6) { continue; } double t = mat[k][k] / mat[i][k]; for (int j = k; j < n + 1; j++) { mat[i][j] = mat[i][j] * t - mat[k][j]; } } }
for (int i = n - 1; i >= 0; i--) { double left = 0; for (int j = i + 1; j < n; j++) { left += ans[j] * mat[i][j]; } if (fabs(mat[i][i]) <= 1e-6) { puts("No Solution"); return 0; } ans[i] = (mat[i][n] - left) / mat[i][i]; }
for (int i = 0; i < n; i++) { printf("%.2lf\n", ans[i]); } return 0; }
|
话说 注意精度问题,判零使用 。