【线性代数入门 | 那忘算8】洛谷P3389 高斯消元(内附行列式教学)
想了想还是单开了一篇,数学王子值得!
专栏指路:《再来一遍一定记住的算法(那些你可能忘记了的算法)》
前置知识:
矩阵:数的集合,一般是方程的系数。
题面:
洛谷P3389 【模板】高斯消元法
说了那么多,其实就是想你解方程组。
1.模拟流程
假设我们有方程:
step 1:写出增广矩阵
(就是把所有数字按照位置写下来啦!)
step 2:定义 r 为将要处理的行,遍历每一列,
再遍历 r 以下的行,找到第一个非零元素
代码这么写:
int r = 1;for (int c = 1;c <= n; c++) { for (int i = r + 1; i <= n; i++) {if ( fabs(a[i][c]) > eps ) {
(fabs 就是专属浮点数的绝对值,eps 是浮点数比较的精度阈值,一般为极小值,eps 约等于 0)
接下来我们要做的就是用 a[i] 行乘以一个倍数,去和 a[r] 行相减,相减完 a[r][c] 就为 0。
double bs = a[r][c] / a[i][c];for (int j = 1; j <= n + 1; j++) {a[r][j] = a[r][j] - a[i][j] * bs;}
比如这个:
第一次循环的 a[i][c] = -3,a[r][c] = 2,bs = 。
整行相减完就变成这样:
然后交换 i 行和 r 行:
swap(a[r], a[i]);
继续寻找 c 列下一个不为 0 的数。
重复以上操作,直到 c 列只有 r 行是不为 0 的。
然后 r++,c++,再次循环。
(顺便一提,每一行第一个不为 0 的数也叫主元,其他的叫自由元,
所以最后的 r - 1 也是主元的个数)
这段的代码:
int r = 1;for (int c = 1;c <= n; c++) { for (int i = r + 1; i <= n; i++) {while ( fabs(a[i][c]) > eps ) { double bs = a[r][c] / a[i][c];for (int j = 1; j <= n + 1; j++) {a[r][j] = a[r][j] - a[i][j] * bs;}swap(a[r], a[i]);}}if (fabs(a[r][c]) > eps) {r++;}}
我懒得全部写了,反正最后消完是这样的:
呵呵,最下面那行可以直接求出一个变量的值,再反带回去消元就好。
代码:
for (int i = n; i >= 1; i--) {for (int j = i + 1; j <= n; j++) {a[i][n + 1] -= x[j] * a[i][j];}x[i] = a[i][n + 1] / a[i][i];if ( fabs(x[i]) < eps ) {x[i] = fabs(x[i]); // 处理浮点误差:将接近 0的值设为 0(避免输出 -0.00)}}
2.洛谷 P3389 整体代码
难点都讲完了,还有个判多解和无解的,我写代码注释里了:
无解:
消完长这样,有行系数全为 0,但最后的常数不为 0。
多解:
消完长这样,有行数全为 0。
就相当于那行不存在了,整个就是个不定方程组。
代码:
#include<bits/stdc++.h>
using namespace std;const double eps = 1e-8;
int n;
double a[110][110], x[110];void gauss() {int r = 1;for (int c = 1;c <= n; c++) { for (int i = r + 1; i <= n; i++) {while ( fabs(a[i][c]) > eps ) { double bs = a[r][c] / a[i][c];for (int j = 1; j <= n + 1; j++) {a[r][j] = a[r][j] - a[i][j] * bs;}swap(a[r], a[i]);}}if (fabs(a[r][c]) > eps) {r++;}}//此时找到了 r-1 个主元 //无解 for (int i = r; i <= n; i++) if( fabs(a[i][n + 1]) > eps ) { cout << "No Solution" << "\n";return ;}//多解 if (r <= n) {cout << "No Solution" << "\n";return ;}for (int i = n; i >= 1; i--) {for (int j = i + 1; j <= n; j++) {a[i][n + 1] -= x[j] * a[i][j];}x[i] = a[i][n + 1] / a[i][i];if ( fabs(x[i]) < eps ) {x[i] = fabs(x[i]); // 处理浮点误差:将接近 0的值设为 0(避免输出 -0.00)}}for (int i = 1; i <= n; i++) {cout << fixed << setprecision(2) << x[i] << "\n";}
}int main() {ios::sync_with_stdio(false);cin.tie(0);cin >> n;for (int i = 1; i <= n; i++) {for (int j = 1; j <= n + 1; j++) {cin >> a[i][j];}}gauss();return 0;
}
3.行列式
如果你只是想知道怎么解方程,那你可以退出了。
现在我们来讲讲高斯消元的其他妙用——
求行列式。
前置知识:
逆序对定义:自行百度
知道求和符号 是啥
行列式:矩阵的值。例:矩阵 的行列式写作
或
。
像这个矩阵:
它的行列式 的公式长这样:
别慌,我来解释下:
其中 是将序列
的元素交换
次得到的一个序列。
(由于每次交换都会产生一个逆序对,所以也可以理解为 是序列
的逆序对数)
举个例子,这是一个 3*3 的矩阵:
我们想算它的行列式,那首先要写出 的全排列:
:逆序对数
,前面的系数为
:逆序对数
,前面的系数为
:逆序对数
,前面的系数为
:逆序对数
,前面的系数为
:逆序对数
,前面的系数为
:逆序对数
,前面的系数为
把系数代入 :
这就是三阶行列式。
还有,有些教程会把行列式写成这样:
其中 是n个元素的对称群(所有排列的集合),
是群里的元素(单个排列),
是排列
的系数(偶排列为
,奇排列为
,奇偶排列定义和前文逆序对个数一样)。
讲了这么多,高斯消元到底能起到什么作用呢??
4.理论支持
前面说过:
矩阵:数的集合,一般是方程的系数。
所以方程的一些操作放在矩阵上也是合理的,just like:
方程组操作 | 矩阵操作 | 目的 |
方程交换位置 | 行交换 | 简化消元顺序 |
方程乘以非零常数 | 行缩放 | 归一化主元(把主元变成 1) |
方程相减消元 | 行替换 | 化为行阶梯型/简化行阶梯型 |
所以之前的变换套在矩阵上也是合理的,高斯消元后的矩阵变成了一个上三角矩阵。
最重要的一点:(上)三角矩阵的行列式等于主对角线元素的乘积!
证明下:
行列式:
对于上三角矩阵,若存在 (位于主对角线下方),则
。
要使乘积 ,必须满足
对所有
成立。
由于 是排列,唯一满足条件的排列是恒等排列(即
)。
因此,行列式中唯一非零的项,是主对角线元素的乘积 ,且其系数为
。
所以,高斯消元能快速求矩阵的行列式。
EXTRA.大整数代码实现
放一个大整数的高斯消元求行列式,想要浮点的按照上面求解方程组的改下就好。
还有些细节写代码注释里了:
void gauss() {//主要做的事情就是辗转相减两行,直到该列除特殊位外都为 0 int r = 1; ans = 1; //r记录特殊位的行(一般情况下 r=c for (int c = 1; c <= n; c++) { //枚举列 for (int i = r + 1; i <= n; i++) { //特殊位下面的行 while (K[i][c]) { //这里可不能改 if __int128 bs = K[r][c] / K[i][c]; //用i行的消r行 for (int j = 1; j <= n; j++) K[r][j] -= K[i][j] * bs;swap(K[r], K[i]); //交换 ans *= -1;//不严格证明下://设 i行值为 a,r行值为 b//第一步操作:i:a+b r:b//第二步操作:i:a+b r:b-(a+b)//第三步操作:i:a+b+(-a) r:-a//将 r的 -a的符号提出来,整个行列式就乘了 -1 }}if (K[r][c] != 0) {r++; //判断多解用的,一般都执行 }}for (int i = 1; i <= n; i++) ans *= K[i][i];//高斯消元完整个矩阵就变成上三角矩阵(主对角行列式),直接乘对角即可 qwrite(ans);
}