当前位置: 首页 > java >正文

【线性代数入门 | 那忘算8】洛谷P3389 高斯消元(内附行列式教学)

想了想还是单开了一篇,数学王子值得!

专栏指路:《再来一遍一定记住的算法(那些你可能忘记了的算法)》

前置知识:

矩阵:数的集合,一般是方程的系数。

题面:

洛谷P3389 【模板】高斯消元法

说了那么多,其实就是想你解方程组。

1.模拟流程

假设我们有方程:

\left\{\begin{matrix} 2x+y-z=8\\ -3x-y+2z=-11\\ -2x+y+2z=-3 \end{matrix}\right.

step 1:写出增广矩阵

\begin{pmatrix} 2 & 1 & -1 & 8\\ -3 & -1 & 2 & -11\\ -2 & 1 & 2 & -3 \end{pmatrix}

(就是把所有数字按照位置写下来啦!)

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;}

比如这个:\begin{pmatrix} 2 & 1 & -1 & 8\\ -3 & -1 & 2 & -11\\ -2 & 1 & 2 & -3 \end{pmatrix}

第一次循环的 a[i][c] = -3,a[r][c] = 2,bs = -\frac{2}{3}

整行相减完就变成这样:

\begin{pmatrix} 0 & \frac{1}{3} & -\frac{7}{3} & \frac{2}{3}\\ -3 & -1 & 2 & -11\\ -2 & 1 & 2 & -3 \end{pmatrix}

然后交换 i 行和 r 行:

swap(a[r], a[i]);

\begin{pmatrix} -3 & -1 & 2 & -11\\0 & \frac{1}{3} & -\frac{7}{3} & \frac{2}{3} \\ -2 & 1 & 2 & -3 \end{pmatrix}

继续寻找 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++;}}

我懒得全部写了,反正最后消完是这样的:

\begin{pmatrix} 1 & 0.5 & -0.5 & 4\\ 0 & 1 & 0.5 & 2.5\\ 0 & 0 & 1 & -1 \end{pmatrix}

呵呵,最下面那行可以直接求出一个变量的值,再反带回去消元就好。

代码:

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 整体代码

难点都讲完了,还有个判多解和无解的,我写代码注释里了:

无解:

\begin{pmatrix} 1 & 0.5 & -0.5 & 4\\ 0 & 1 & 0.5 & 2.5\\ 0 & 0 & 0 & -1 \end{pmatrix}

消完长这样,有行系数全为 0,但最后的常数不为 0。

多解:

\begin{pmatrix} 1 & 0.5 & -0.5 & 4\\ 0 & 1 & 0.5 & 2.5\\ 0 & 0 & 0 & 0 \end{pmatrix}

消完长这样,有行数全为 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.行列式

如果你只是想知道怎么解方程,那你可以退出了。

现在我们来讲讲高斯消元的其他妙用——

求行列式

前置知识:

逆序对定义:自行百度

知道求和符号 \sum​ 是啥

 行列式:矩阵的值。例:矩阵 A​ 的行列式写作 |A|​ 或 det(A)​。

像这个矩阵:

A=\begin{vmatrix} a_{11} & a_{12} & ... & a_{1n}\\ a_{21} & a_{22} & ... & a_{2n}\\ ... & ... & ... & ...\\ a_{n1} & a_{n2} & ... & a_{nn} \end{vmatrix}

它的行列式 D​ 的公式长这样:

D=\sum_{k}^{}(-1)^ka_{1k_1}a_{2k_2}...a_{nk_n}

别慌,我来解释下:

其中 k1,k2,...,kn​ 是将序列 1,2,...,n​ 的元素交换 k​ 次得到的一个序列。

(由于每次交换都会产生一个逆序对,所以也可以理解为 k​ 是序列 k1,k2,...,kn​ 的逆序对数

举个例子,这是一个 3*3 的矩阵:

A=\begin{vmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{vmatrix}

我们想算它的行列式,那首先要写出 1,2,3​ 的全排列:

1,2,3​ :逆序对数 k=0​,前面的系数为 (-1)^k=(-1)^0=1

1,3,2​ :逆序对数 k=1​,前面的系数为 (-1)^k=(-1)^1=-1

2,1,3​ :逆序对数 k=1​,前面的系数为 (-1)^k=(-1)^1=-1

2,3,1​ :逆序对数 k=2​,前面的系数为 (-1)^k=(-1)^2=1

3,1,2​ :逆序对数 k=2​,前面的系数为 (-1)^k=(-1)^2=1

3,2,1​ :逆序对数 k=3​,前面的系数为 (-1)^k=(-1)^3=-1

把系数代入 D=\sum_{k}^{}(-1)^ka_{1k_1}a_{2k_2}...a_{nk_n}​:

D=a_{11}a_{22}a_{33}\,-\,a_{11}a_{23}a_{32}\,-\,a_{12}a_{21}a_{33}\,+\,a_{12}a_{23}a_{31}\,+\,a_{13}a_{21}a_{32}\,-\,a_{13}a_{22}a_{31}

这就是三阶行列式

还有,有些教程会把行列式写成这样:

det(A)=\sum_{\sigma \in S_n}^{}sgn(\sigma )a_{1\sigma (1)}a_{2\sigma (2)}...a_{n\sigma (n)}

其中 S_n 是n个元素的对称群(所有排列的集合),

\sigma 是群里的元素(单个排列),

sign(\sigma ) 是排列 \sigma 的系数(偶排列为 1,奇排列为 -1,奇偶排列定义和前文逆序对个数一样)。

讲了这么多,高斯消元到底能起到什么作用呢??

4.理论支持

前面说过:

矩阵:数的集合,一般是方程的系数。

所以方程的一些操作放在矩阵上也是合理的,just like:

方程组操作矩阵操作目的
方程交换位置行交换简化消元顺序
方程乘以非零常数行缩放归一化主元(把主元变成 1)
方程相减消元行替换化为行阶梯型/简化行阶梯型

所以之前的变换套在矩阵上也是合理的,高斯消元后的矩阵变成了一个上三角矩阵

\begin{pmatrix} 1 & 0.5 & -0.5 & 4\\ 0 & 1 & 0.5 & 2.5\\ 0 & 0 & 1 & -1 \end{pmatrix}

最重要的一点:(上)三角矩阵的行列式等于主对角线元素的乘积!

证明下:

行列式:det(A)=\sum_{\sigma \in S_n}^{}sgn(\sigma )a_{1\sigma (1)}a_{2\sigma (2)}...a_{n\sigma (n)}

对于上三角矩阵,若存在 i>\sigma (i)(位于主对角线下方),则 a_{i\sigma (i)} = 0

要使乘积 a_{1\sigma (1)}a_{2\sigma (2)}...a_{n\sigma (n)},必须满足 \sigma (i)\geq i 对所有 i 成立。

由于 \sigma 是排列,唯一满足条件的排列是恒等排列(即 \sigma (i)= i)。

因此,行列式中唯一非零的项,是主对角线元素的乘积 a_{11}a_{22}...a_{nn},且其系数为 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);
}

http://www.xdnf.cn/news/19083.html

相关文章:

  • web3简介
  • 屏随人动+视觉魔方+多样主题+智能留言,涂鸦Wukong AI 2.0助力打造爆款带屏云台相机
  • DVWA靶场通关笔记-命令执行(Impossible级别)
  • 如何制作手感良好的移动算法?
  • 【视频讲解】R语言海七鳃鳗性别比分析:JAGS贝叶斯分层逻辑回归MCMC采样模型应用
  • GPT-Realtime架构与Token成本控制深度解析
  • 解析DB-GPT项目中三个 get_all_model_instances 方法的区别
  • 考研数据结构Part3——二叉树知识点总结
  • 大数据毕业设计选题推荐:基于北京市医保药品数据分析系统,Hadoop+Spark技术详解
  • useEffect用法
  • 将2D基础模型(如SAM/SAM2)生成的2D语义掩码通过几何一致性约束映射到3D高斯点云
  • 告别K8s部署繁琐!用KubeOperator可视化一键搭建生产级集群
  • 数据结构 02(线性:顺序表)
  • aggregating英文单词学习
  • 数字人 + 矩阵聚合系统源码搭建与定制化开发
  • Python 轻量级 HTML 解析器 - lxml入门教程
  • 通过Kubernetes安装mysql5服务
  • 深入解析Qt节点编辑器框架:数据流转与扩展机制(三)
  • 4. LangChain4j 模型参数配置超详细说明
  • 机器学习回顾——线性回归
  • Redis红锁(RedLock)解密:分布式锁的高可用终极方案
  • DBeaver中禁用PostgreSQL SSL的配置指南
  • 【性能优化】Unity 渲染优化全解析:Draw Call、Batch、SetPass 与批处理技术
  • 【Django】首次创建Django项目初始化
  • “帕萨特B5钳盘式制动器结构设计三维PROE模型7张CAD图纸PDF图“
  • 人工智能基础概念
  • 秋招笔记-8.28
  • 总结:在工作场景中的应用。(Excel)
  • Dify学习
  • 响应式编程框架Reactor【1】