草庐IT

「学习笔记」高斯消元

Keven-He 2023-03-28 原文

「学习笔记」高斯消元

点击查看目录

目录

注意：本文内的高斯消元均是「高斯-约旦消元法（\(\text{Gauss-Jordan Elimination}\)）」。

突然发现有的电脑好像看不了增广矩阵里的竖线，那就凑合看吧（

算法

思路

我们举个例子：

\[\begin{cases} x_2+x_3&=3\\ x_1+2x_3&=9\\ 2x_1+x_2+x_3&=13 \end{cases} \]

我们把它的所有系数化为矩阵 \(A\)：

\[\left[\begin{matrix} 0&1&1\\ 1&0&2\\ 2&1&1\\ \end{matrix} \middle| \begin{matrix} 3 \\ 9 \\ 13 \end{matrix} \right] \]

然后开始推：

  1. 首先遍历第 \(i\) 列的第 \(i\) 行至 \(n\) 行（因为前 \(i-1\) 行已经形成阵型了，不能破坏），找出这一列最大的系数，并将那一行交换至第 \(i\) 行。（当前 \(i=1\)）

\[\left[\begin{matrix} 2&1&1\\ 1&0&2\\ 0&1&1\\ \end{matrix} \middle| \begin{matrix} 13 \\ 9 \\ 3 \\ \end{matrix} \right] \]

  1. 然后判断解，即判断 \(A_{i,i}\) 是否为零。若为零则无唯一解，否则有唯一解。（这个例子显然有唯一解）

  2. 系数化为一，即将 \(A_{i,j}(i<j\le n+1)\) 除以 \(A_{i,i}\)，再将 \(A_{i,i}\) 除以自己得 \(1\)。

\[\left[\begin{matrix} 1&0.5&0.5\\ 1&0&2\\ 0&1&1\\ \end{matrix} \middle| \begin{matrix} 6.5 \\ 9 \\ 3 \\ \end{matrix} \right] \]

  1. 然后用加减消元法消掉一部分元：

\[\left[\begin{matrix} 1&0.5&0.5\\ 0&-0.5&1.5\\ 0&1&1\\ \end{matrix} \middle| \begin{matrix} 6.5 \\ 2.5 \\ 3 \\ \end{matrix} \right] \]

  1. 重复上述步骤：

\[\left[\begin{matrix} 1&0.5&0.5\\ 0&-0.5&1.5\\ 0&1&1\\ \end{matrix}\middle|\begin{matrix} 6.5 \\ 2.5 \\ 3 \\ \end{matrix}\right] \Longrightarrow \left[\begin{matrix} 1&0.5&0.5\\ 0&1&1\\ 0&-0.5&1.5\\ \end{matrix}\middle|\begin{matrix} 6.5 \\ 3 \\ 2.5 \\ \end{matrix}\right] \Longrightarrow \left[\begin{matrix} 1&0&0\\ 0&1&1\\ 0&0&2\\ \end{matrix}\middle|\begin{matrix} 5 \\ 3 \\ 4 \\ \end{matrix}\right] \Longrightarrow \left[\begin{matrix} 1&0&0\\ 0&1&1\\ 0&0&1\\ \end{matrix}\middle|\begin{matrix} 5 \\ 3 \\ 2 \\ \end{matrix}\right] \Longrightarrow \left[\begin{matrix} 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{matrix}\middle|\begin{matrix} 5 \\ 1 \\ 2 \\ \end{matrix}\right] \]

所有操作步骤完成，最后得到解：

\[\begin{cases} x_1=5\\ x_2=1\\ x_3=2\\ \end{cases} \]

时间复杂度 \(\Theta(n^3)\)。

注意：输入不要用快读，虽然有的题输入只有整数，但很多题都不是！

代码

点击查看代码 
inline void Gauss(){
	_for(i,1,n){
		ll mx=i;
		_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
		_for(j,1,n+1)swap(a[i][j],a[mx][j]);
		if(fabs(a[i][i])<eps){puts("No Solution");return;}
		_for(j,i+1,n+1)a[i][j]/=a[i][i];
		a[i][i]=1;
		_for(j,1,n){
			if(i==j)continue;
			_for(k,i+1,n+1)a[j][k]-=a[j][i]*a[i][k];
		}
	}
	_for(i,1,n)printf("%.2lf\n",a[i][n+1]);
}

例题:[SDOI2006]线性方程组

思路

本题主要难点在于如何区分无解和无限组解。

举个例子：

\[A= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&2&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 5 \\ \end{matrix}\right], B= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&2&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 7 \\ \end{matrix}\right], C= \left[\begin{matrix} 1&0&1\\ 0&2&1\\ 1&1&2\\ \end{matrix}\middle|\begin{matrix} 2 \\ 3 \\ 5 \\ \end{matrix}\right] \]

  • 在 \(A\) 中，第三行的前 \(3\) 项是前两行相加的结果，第 \(4\) 项也是，所以有用的行只有两个（这个数量称为秩），有一个数是不确定的，所以有不唯一解。

  • 在 \(B\) 中，第三行的前 \(3\) 项是前两行相加的结果，但第 \(4\) 项不是，会出现 \(0=2\) 的情况，所以无解。

  • 在 \(C\) 中，它的秩是 \(3\)，所以有唯一解：

\[\begin{cases} x_1=-1\\ x_2=0\\ x_3=3\\ \end{cases} \]

归纳一下：

  • 若秩不等于未知数的数量且消元后出现 \(\sum_{j=1}^{n}a_{i,j}=0\) 且 \(a_{i,n+1}\neq0\) 的情况，则无解。

  • 若秩不等于未知数的数量且消元后出现 \(\sum_{j=1}^{n+1}a_{i,j}=0\) 的情况，则有无数组解。

  • 否则有唯一解。

代码

点击查看代码 
const ll N=110,inf=1ll<<40;
const double eps=1e-9;
ll n,l=1;db a[N][N];
namespace SOLVE{
	inline void Gauss(){
		_for(i,1,n){
			ll mx=l;
			_for(j,l+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
			if(fabs(a[mx][i])<eps)continue;
			swap(a[mx],a[l]);
			_for(j,i+1,n+1)a[l][j]/=a[l][i];
			a[l][i]=1;
			_for(j,1,n){
				if(j==l)continue;
				_for(k,i+1,n+1)a[j][k]-=a[l][k]*a[j][i];
			}
			++l;
		}
	}
	inline void Jie(){
		if(l<=n){
			_for(i,l,n)if(fabs(a[i][n+1])>eps){puts("-1");return;}
			puts("0");return;
		}
		_for(i,1,n){
			printf("x%lld=%.2lf\n",i,a[i][n+1]);
		}
	}
	inline void In(){
		scanf("%lld",&n);
		_for(i,1,n)_for(j,1,n+1)scanf("%lf",&a[i][j]);
		Gauss();
		Jie();
		return;
	}
}

练习题

[JSOI2008]球形空间产生器

思路

设球心为 \(p(p_1,p_2,\cdots,p_n)\)，可列式子：

\[\begin{aligned} r=&(a_{1,1}-p_1)^2+(a_{1,2}-p_2)^2+\cdots+(a_{1,n}-p_n)^2\\ =&(a_{2,1}-p_1)^2+(a_{2,2}-p_2)^2+\cdots+(a_{2,n}-p_n)^2\\ =&\cdots\\ =&(a_{n+1,1}-p_1)^2+(a_{n+1,2}-p_2)^2+\cdots+(a_{n+1,n}-p_n)^2\\ \end{aligned} \]

拆一下得：

\[\begin{aligned} r=&(a_{1,1}^2-2a_{1,1}p_1+p_1^2)+(a_{1,2}^2-2a_{1,2}p_2+p_2^2)+\cdots+(a_{1,n}^2-2a_{1,n}p_n+p_n^2)\\ =&(a_{2,1}^2-2a_{2,1}p_1+p_1^2)+(a_{2,2}^2-2a_{2,2}p_2+p_2^2)+\cdots+(a_{2,n}^2-2a_{2,n}p_n+p_n^2)\\ =&\cdots\\ =&(a_{n+1,1}^2-2a_{n+1,1}p_1+p_1^2)+(a_{n+1,2}^2-2a_{n+1,2}p_2+p_2^2)+\cdots+(a_{n+1,n}^2-2a_{n+1,n}p_n+p_n^2)\\ \end{aligned} \]

化简得：

\[\begin{aligned} r=&\sum_{i=1}^{n}a_{1,i}^2-\sum_{i=1}^{n}2a_{1,i}p_i\\ =&\sum_{i=1}^{n}a_{2,i}^2-\sum_{i=1}^{n}2a_{2,i}p_i\\ =&\cdots\\ =&\sum_{i=1}^{n}a_{n+1,i}^2-\sum_{i=1}^{n}2a_{n+1,i}p_i\\ \end{aligned} \]

即：

\[\begin{aligned} \sum_{i=1}^{n}a_{1,i}^2=&r+\sum_{i=1}^{n}2a_{1,i}p_i\\ \sum_{i=1}^{n}a_{2,i}^2=&r+\sum_{i=1}^{n}2a_{2,i}p_i\\ \cdots=&r+\cdots\\ \sum_{i=1}^{n}a_{n+1,i}^2=&r+\sum_{i=1}^{n}2a_{n+1,i}p_i\\ \end{aligned} \]

发现就是个方程，但是这个 \(r\) 混在里面很难受。那么我们把它当成一个系数为 \(1\) 的未知数，则原式可以化成如下增广矩阵：

\[\left[ \begin{matrix} 2a_{1,1}&2a_{1,2}&\cdots&2a_{1,n}&1\\\\ 2a_{2,1}&2a_{2,2}&\cdots&2a_{2,n}&1\\\\ \cdots&\cdots&\cdots&\cdots&1\\\\ 2a_{n+1,1}&2a_{n+1,2}&\cdots&2a_{n+1,n}&1 \end{matrix} \middle| \begin{matrix} \sum_{i=1}^{n}a_{1,i}^2\\\\ \sum_{i=1}^{n}a_{2,i}^2\\\\ \cdots\\\\ \sum_{i=1}^{n}a_{n+1,i}^2\\ \end{matrix} \right] \]

就可以直接高斯消元了。

Latex 当算草纸真好用。

代码

点击查看代码 
const ll N=20,inf=1ll<<40;
const double eps=1e-6;
ll n;ldb a[N][N];
namespace SOLVE{
	inline ll rnt(){
		ll x=0,w=1;char c=getchar();
		while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
		while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
		return x*w;
	}
	inline void Gauss(){
		_for(i,1,n){
			ll mx=i;
			_for(j,i+1,n)if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
			swap(a[mx],a[i]);
			_for(j,i+1,n+1)a[i][j]/=a[i][i];
			a[i][i]=1.0;
			_for(j,1,n){
				if(i==j)continue;
				for_(k,n+1,i)a[j][k]-=a[j][i]*a[i][k];
			}
		}
	}
	inline void In(){
		scanf("%lld",&n),++n;
		_for(i,1,n){
			_for(j,1,n-1){
				scanf("%Lf",&a[i][j]);
				a[i][n+1]+=a[i][j]*a[i][j];
				a[i][j]*=2.0;
			}
			a[i][n]=1;
		}
		Gauss();
		_for(i,1,n-1)printf("%.3Lf ",a[i][n+1]);
		return;
	}
}

[USACO10HOL]Driving Out the Piggies G

思路

好题，还对高斯消元能出什么题没什么认知就 刷新了我对高斯消元能出什么题的认知。

我们设到节点 \(i\) 的期望步数为 \(x_i\)，节点 \(i\) 在图上的出度为 \(d_i\)。

那么显然：

\[x_i=[i=1]+\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times \frac{1}{d_j} \]

显然 \(x_i\times\frac{p}{q}\) 为实际爆炸的概率，那么考虑如何求 \(x_i\)。我们进行一个移项：

\[x_i-\sum_{(i,j)\in\mathbb{E}}x_j\times(1-\frac{p}{q})\times \frac{1}{d_j}=[i=1] \]

直接冲一个高斯消元即可。

代码

点击查看代码 
const ll N=310,inf=1ll<<40;
const double eps=1e-13;
ll n,m,tu[N][N];ldb p,q,d[N],a[N][N];
namespace SOLVE{
	inline ll rnt(){
		ll x=0,w=1;char c=getchar();
		while(!isdigit(c)){if(c=='-')w=-1;c=getchar();}
		while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
		return x*w;
	}
	inline void Pre(){
		a[1][n+1]=1;
		_for(i,1,n){
			_for(j,1,n){
				if(i==j)a[i][j]=1.0;
				else if(tu[i][j])a[i][j]=-(1.0-p/q)/d[j];
			}
		}
	}
	inline void Gauss(){
		_for(i,1,n){
			ll mx=i;
			_for(j,i+1,n)if(a[mx][i]<a[j][i])mx=j;
			swap(a[mx],a[i]);
			_for(j,i+1,n+1)a[i][j]/=a[i][i];
			a[i][i]=1;
			_for(j,1,n){
				if(i==j)continue;
				for_(k,n+1,i)a[j][k]-=a[i][k]*a[j][i];
			}
		}
	}
	inline void In(){
		scanf("%lld%lld%Lf%Lf",&n,&m,&p,&q);
		_for(i,1,m){
			ll x=rnt(),y=rnt();
			d[x]+=1.0,d[y]+=1.0;
			tu[x][y]=tu[y][x]=1;
		}
		Pre(),Gauss();
		_for(i,1,n)printf("%.9Lf\n",a[i][n+1]*p/q);
		return;
	}
}

\[\Huge{\mathfrak{The\ End}} \]

笔记学习ampspanclass算法与数据结构

有关「学习笔记」高斯消元的更多相关文章

  1. LC滤波器设计学习笔记（一）滤波电路入门 - 2

    目录前言滤波电路科普主要分类实际情况单位的概念常用评价参数函数型滤波器简单分析滤波电路构成低通滤波器RC低通滤波器RL低通滤波器高通滤波器RC高通滤波器RL高通滤波器部分摘自《LC滤波器设计与制作》，侵权删。前言最近需要学习放大电路和滤波电路，但是由于只在之前做音乐频谱分析仪的时候简单了解过一点点运放，所以也是相当从零开始学习了。滤波电路科普主要分类滤波器：主要是从不同频率的成分中提取出特定频率的信号。有源滤波器：由RC元件与运算放大器组成的滤波器。可滤除某一次或多次谐波，最普通易于采用的无源滤波器结构是将电感与电容串联，可对主要次谐波（3、5、7）构成低阻抗旁路。无源滤波器：无源滤波器，又称

  2. CAN协议的学习与理解 - 2

    最近在学习CAN，记录一下，也供大家参考交流。推荐几个我觉得很好的CAN学习，本文也是在看了他们的好文之后做的笔记首先是瑞萨的CAN入门，真的通透；秀！靠这篇我竟然2天理解了CAN协议！实战STM32F4CAN！原文链接：https://blog.csdn.net/XiaoXiaoPengBo/article/details/116206252CAN详解（小白教程）原文链接：https://blog.csdn.net/xwwwj/article/details/105372234一篇易懂的CAN通讯协议指南1一篇易懂的CAN通讯协议指南1-知乎(zhihu.com)视频推荐CAN总线个人知识总

  3. 深度学习部署：Windows安装pycocotools报错解决方法 - 2

    深度学习部署：Windows安装pycocotools报错解决方法1.pycocotools库的简介2.pycocotools安装的坑3.解决办法更多Ai资讯：公主号AiCharm本系列是作者在跑一些深度学习实例时，遇到的各种各样的问题及解决办法，希望能够帮助到大家。ERROR:Commanderroredoutwithexitstatus1:'D:\Anaconda3\python.exe'-u-c'importsys,setuptools,tokenize;sys.argv[0]='"'"'C:\\Users\\46653\\AppData\\Local\\Temp\\pip-instal

  4. ruby - 我正在学习编程并选择了 Ruby。我应该升级到 Ruby 1.9 吗？ - 2

    我完全不是程序员，正在学习使用Ruby和Rails框架进行编程。我目前正在使用Ruby1.8.7和Rails3.0.3，但我想知道我是否应该升级到Ruby1.9，因为我真的没有任何升级的“遗留”成本。缺点是什么？我是否会遇到与普通gem的兼容性问题，或者甚至其他我不太了解甚至无法预料的问题？ 最佳答案 你应该升级。不要坚持从1.8.7开始。如果您发现不支持1.9.2的gem，请避免使用它们(因为它们很可能不被维护)。如果您对gem是否兼容1.9.2有任何疑问，您可以在以下位置查看:http://www.railsplugins.or

  5. ruby - 我如何学习 ruby​​ 的正则表达式？ - 2

    如何学习ruby​​的正则表达式？(对于假人) 最佳答案 http://www.rubular.com/在Ruby中使用正则表达式时是一个很棒的工具，因为它可以立即将结果可视化。 关于ruby-我如何学习ruby​​的正则表达式？，我们在StackOverflow上找到一个类似的问题： https://stackoverflow.com/questions/1881231/

  6. 深度学习12. CNN经典网络 VGG16 - 2

    深度学习12.CNN经典网络VGG16一、简介1.VGG来源2.VGG分类3.不同模型的参数数量4.3x3卷积核的好处5.关于学习率调度6.批归一化二、VGG16层分析1.层划分2.参数展开过程图解3.参数传递示例4.VGG16各层参数数量三、代码分析1.VGG16模型定义2.训练3.测试一、简介1.VGG来源VGG（VisualGeometryGroup）是一个视觉几何组在2014年提出的深度卷积神经网络架构。VGG在2014年ImageNet图像分类竞赛亚军，定位竞赛冠军；VGG网络采用连续的小卷积核（3x3）和池化层构建深度神经网络，网络深度可以达到16层或19层，其中VGG16和VGG

  7. 机器学习——时间序列ARIMA模型(四)：自相关函数ACF和偏自相关函数PACF用于判断ARIMA模型中p、q参数取值 - 2

    文章目录1、自相关函数ACF2、偏自相关函数PACF3、ARIMA(p,d,q)的阶数判断4、代码实现1、引入所需依赖2、数据读取与处理3、一阶差分与绘图4、ACF5、PACF1、自相关函数ACF自相关函数反映了同一序列在不同时序的取值之间的相关性。公式：ACF(k)=ρk=Cov(yt,yt−k)Var(yt)ACF(k)=\rho_{k}=\frac{Cov(y_{t},y_{t-k})}{Var(y_{t})}ACF(k)=ρk​=Var(yt​)Cov(yt​,yt−k​)​其中分子用于求协方差矩阵，分母用于计算样本方差。求出的ACF值为[-1,1]。但对于一个平稳的AR模型，求出其滞

  8. Unity Shader 学习笔记（5）Shader变体、Shader属性定义技巧、自定义材质面板 - 2

    写在之前Shader变体、Shader属性定义技巧、自定义材质面板，这三个知识点任何一个单拿出来都是一套知识体系，不能一概而论，本文章目的在于将学习和实际工作中遇见的问题进行总结，类似于网络笔记之用，方便后续回顾查看，如有以偏概全、不祥不尽之处，还望海涵。1、Shader变体先看一段代码......Properties{ [KeywordEnum(on,off)]USL_USE_COL("IsUseColorMixTex?",int)=0 [Toggle(IS_RED_ON)]_IsRed("IsRed?",int)=0}......//中间省略，后续会有完整代码 #pragmamulti_c

  9. Tcl脚本入门笔记详解（一） - 2

    TCL脚本语言简介•TCL（ToolCommandLanguage）是一种解释执行的脚本语言（ScriptingLanguage）,它提供了通用的编程能力：支持变量、过程和控制结构；同时TCL还拥有一个功能强大的固有的核心命令集。TCL经常被用于快速原型开发，脚本编程，GUI和测试等方面。•实际上包含了两个部分：一个语言和一个库。首先，Tcl是一种简单的脚本语言，主要使用于发布命令给一些互交程序如文本编辑器、调试器和shell。由于TCL的解释器是用C\C++语言的过程库实现的，因此在某种意义上我们又可以把TCL看作C库，这个库中有丰富的用于扩展TCL命令的C\C++过程和函数，所以，Tcl是

  10. ruby-on-rails - 这个 C 和 PHP 程序员如何学习 Ruby 和 Rails？ - 2

    按照目前的情况，这个问题不适合我们的问答形式。我们希望答案得到事实、引用或专业知识的支持，但这个问题可能会引发辩论、争论、投票或扩展讨论。如果您觉得这个问题可以改进并可能重新打开，visitthehelpcenter指导。关闭9年前。我来自C、php和bash背景，很容易学习，因为它们都有相同的C结构，我可以将其与我已经知道的联系起来。然后2年前我学了Python并且学得很好，Python对我来说比Ruby更容易学。然后从去年开始，我一直在尝试学习Ruby，然后是Rails，我承认，直到现在我还是学不会，讽刺的是那些打着简单易学的烙印，但是对于我这样一个老练的程序员来说，我只是无法将它

随机推荐