浅谈高斯消元的实现和简单应用
一、高斯消元的原理
对于n元的m个线性方程组成的方程组,我们将其以矩阵的形式记录下来:
a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn
然后进行初等行列变换,尝试构造出一个上三角矩阵,逐步使系数不为零的项减少;
等最后只剩下一个系数不为零时,进行回代,逐步求出已知解。(详解过程咨询小学老师)
二、高斯消元的实现
老老实实的回代代码参见其他人的博客,这里介绍一种比较毒瘤的不回代暴力消元法:
1.Process
对于每个方程,按照一定的规则(后话)挑选一个主元,记录该主元对应第几个方程,然后用初等行列变换消去其他所有该元的系数;
最后我们得到的是一个每行只有一个数不为零,每列只有一个数不为零的鬼畜矩阵(自己脑补)
此时令ans向量对应的数字出去该行的非零系数,即为对应该元的解。
2.Code
设a数组为已经记录系数的数组(格式见上方),那么a应该是n行n+1列的,最后一列代表方程的常数项;
设w数组记录每个方程的主元是第几项,v数组记录答案;
当多解时输出“Multiple solutions”,无解时输出”No solution”;
double a[max_n][max_n+1],v[max_n];int w[max_n]; void gauss(){double eps=1e-6;for(int i=1;i<=n;++i){ //Enumerate the equation;int p=0; //Record the position of the largest number; double mx=0; //Recording the largest number;for(int j=1;j<=n;++j)if(fabs(a[i][j])-eps>mx){mx=fabs(a[i][j]);p=j; //fabs() returns the absolute value of float; }if(!p){if(fabs(a[i][n+1]<eps))printf("Multiple solutions");else printf("No solution");return; }w[i]=p;for(int j=1;j<=n;++j)if(i!=j){ //other equationsdouble t=a[j][p]/a[i][p];for(int k=1;k<=n+1;++k) //n+1 is importanta[j][k]-=a[i][k]*t;}}for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
}
3.notice
(1)精度的设置
众所周知浮点数是有精度丢失的,在高斯消元中,精度的选择要依题目而定,精度过低会导致系数较小的数被误认为系数为零,而精度过高也有可能会导致误差大于精度,导致本应该系数消为0的项误认为系数不为零,所以精度的选择是很哲学的问题,要注意。
推荐范围:1e-4到1e-10
(2)主元的选取原则
接着(1)说,我们知道,用浮点数是有精度丢失的,那么让一个较大的数除以一个极小的数误差自然大的可想而知,所以我们想得到在精度允许的条件下系数最大的主元,所以对于每个方程,我们都应该选择最大系数的元作为主元。
(3)在模2意义下的高斯消元
使用bitset优化运行时间,详见相关应用中第三个例题的讲解;
三、相关应用
这里给出高斯消元的几道基础题目,难度适合初学者。
1.[Luogu P3389]【模板】高斯消元
Description
给定一个线性方程组,对其求解
输入格式:
第一行,一个正整数 n
第二至 n+1行,每行 n+1个整数,为 a1,a2⋯an和 b,代表一组方程。
输出格式:
共n行,每行一个数,第 i行为 xi(保留2位小数)
如果不存在唯一解,在第一行输出"No Solution".
Solution
如上所述。
Code
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;const int max_n=110;
double a[max_n][max_n+1],v[max_n];
int n,w[max_n]; inline int rd(){int x=0;bool f=0;char c=getchar();while(!isdigit(c)){if(c=='-') f=1;c=getchar();}while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}return f?-x:x;
}void gauss(){double eps=1e-6;for(int i=1;i<=n;++i){//enumerate the equation;int p=0; //Record the position of the largest number; double mx=0; //Recording the largest number;for(int j=1;j<=n;++j)if(fabs(a[i][j])-eps>mx){mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; }if(!p){printf("No Solution");return; }w[i]=p;for(int j=1;j<=n;++j)if(i!=j){ //other equationsdouble t=a[j][p]/a[i][p];for(int k=1;k<=n+1;++k)//n+1 is importanta[j][k]-=a[i][k]*t;}}for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];for(int i=1;i<=n;++i) printf("%.2lf\n",v[i]);
}int main(){n=rd();for(int i=1;i<=n;++i)for(int j=1;j<=n+1;++j)a[i][j]=rd();gauss();return 0;
}
2.[BZOJ 1013][JSOI 2008] 球形空间产生器sphere
详解参考我的随笔:http://www.cnblogs.com/COLIN-LIGHTNING/p/8982341.html
转载于:https://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html
浅谈高斯消元的实现和简单应用相关推荐
- 解线性方程组——高斯消元の板子
ATP记得它在很久以前看过一点点高斯消元的东西然后做过一点点题目..但是当时实在是太zz了所以本来就没有很懂这个东西现在更是忘得差不多了.. 所以现在就当重新学一遍了QwQ 一点口胡的解释 高斯消元. ...
- Rocksdb Ribbon Filter : 结合 XOR-filter 以及 高斯消元算法 实现的 高效filter
文章目录 前言 XOR-filter 实现原理 xor filter 的构造原理 xor filter 构造总结 XOR-filter 和 ADD-filter对比 XOR-filter 在计算上的优 ...
- poj 1681 Painter#39;s Problem(高斯消元)
http://poj.org/problem? id=1681 求最少经过的步数使得输入的矩阵全变为y. 思路:高斯消元求出自由变元.然后枚举自由变元,求出最优值. 注意依据自由变元求其它解及求最优值 ...
- AC自动机 + 概率dp + 高斯消元 --- HDU 5955 or 2016年沈阳icpc H [AC自动机 + 概率dp + 高斯消元]详解
题目链接 题目大意: 就是有NNN个人,每个人都会猜一个长度为LLL的只包含{1,2,3,4,5,6}\{1,2,3,4,5,6\}{1,2,3,4,5,6}的序列,现在裁判开始投掷骰子,并且把每次的 ...
- ICPC 2005 hangzhou Generator (UVA1358)KMP + 期望DP / 高斯消元
整理的算法模板合集: ACM模板 点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 Generator Weblink https://www.luogu.com.cn/problem/ ...
- 2020 ACM / ICPC 济南 A Matrix Equation (高斯消元、乘法原理)
整理的算法模板合集: ACM模板 点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 题目链接 给你定义两种 010101 矩阵上的运算: Xi,j×Yi,j=(∑k=1NXi,kYk,j ...
- luogu P4035 [JSOI2008]球形空间产生器(高斯消元 / 模拟退火)
整理的算法模板合集: ACM模板 点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 数据范围只开到了10,而且是经典的力学结构,所以我们可以用模拟退火,可以做一下 nnn 维的正交分解h ...
- BZOJ 2707: [SDOI2012]走迷宫 [高斯消元 scc缩点]
2707: [SDOI2012]走迷宫 题意:求s走到t期望步数,\(n \le 10^4\),保证\(|SCC| \le 100\) 求scc缩点,每个scc高斯消元,scc之间直接DP 注意每次清 ...
- [BZOJ 3143][Hnoi2013]游走(高斯消元+期望)
Description 一个无向连通图,顶点从1编号到N,边从1编号到M. 小Z在该图上进行随机游走,初始时小Z在1号顶点,每一步小Z以相等的概率随机选 择当前顶点的某条边,沿着这条边走到下一个顶点, ...
最新文章
- Android 逐帧动画(Frame)
- 转型AI成功几率有几分?太真实了......
- 阿里巴巴副总裁陈丽娟:我对阿里云产品生态的思考 | 云原生加速器观点
- Unable to locate package php5-curl
- boost::lexical_cast用法的测试程序
- C#开发笔记之19-如何用C#实现优雅的Json解析(序列化/反序列化)方案?
- 用bool函数判断int类型相加溢出_Go是强类型语言,不支持隐式类型转换,那该怎么办?...
- group by 深入总结
- Oracle傻瓜手册
- IDEA 创建java项目
- 软件开发中JAVA编程语言的应用
- 华为确认53岁高管丁耘骤逝:执掌最大营收业务,东南大学毕业,在职已超26年...
- 电脑激活Office时出现异常,激活界面白屏或提示无法与服务器
- Nvidia TX2 串口使用
- 2022年30本新年书单(要么旅行,要么读书,身体和灵魂总有一个在路上)
- View事件分发机制分析
- 【机器学习】异常检测算法之(KNN)-K Nearest Neighbors
- 【190319】VC++ C/S结构视频聊天软件源码源代码
- linux断开网卡,Linux 无线网络断开的解决方案
- Ansoft安装报错reg_ansys.exe ERROR:Error3221227010
热门文章
- Rstudio调用plot()函数时,出现错误的处理方法
- 简单的一个月之设计实现内容1
- 记一次sql优化之索引的引用
- 20180307:python接口测试时json的传参与解析区分
- HTML 中的特殊字符
- select * 和select 所有字段的区别
- [修正] Berlin 10.1 支持 iPhone 4 (iOS v7.x)
- mybatis判断集合为空或者元素个数为零
- 【观点】从曾成杰案看民间金融的高风险与银行缺失的机制创新
- Asp.net 调用mysql存储过程参数传中文乱码!