博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
浅谈高斯消元的实现和简单应用
阅读量:5126 次
发布时间:2019-06-13

本文共 3058 字,大约阅读时间需要 10 分钟。

一、高斯消元的原理

对于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]

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
#include
#include
#include
#include
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 equations double t=a[j][p]/a[i][p]; for(int k=1;k<=n+1;++k)//n+1 is important a[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

详解参考我的随笔:

转载于:https://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html

你可能感兴趣的文章
solr的一些查询语法
查看>>
你所不了解的五条面试忠告
查看>>
每日一题20180330-Linux
查看>>
从零开始学习Hadoop--第2章 第一个MapReduce程序
查看>>
unity3d 捕获系统日志,来处理一些问题
查看>>
Android的Recovery中font_10x10.h字库文件制作
查看>>
SpringBoot+FreeMarker开发word文档下载,预览
查看>>
JVM 内存设置大小(Xms Xmx PermSize MaxPermSize 区别)
查看>>
STL中 map 和 multimap
查看>>
discuz 标签详解
查看>>
微信公众账户模拟登陆后的一系列操作
查看>>
Mac远程连接服务器
查看>>
使用python爬取东方财富网机构调研数据
查看>>
贝叶斯理论基础理解
查看>>
2018java最新面试题
查看>>
PHP编写命令行脚本和后台运行程序的注意事项
查看>>
php 换行 PHP_EOL变量
查看>>
JS中关于clientWidth、offsetWidth、scrollWidth
查看>>
Telerik Reporting之生成报表
查看>>
不要怂!就是干!
查看>>