浅谈高斯消元的实现和简单应用

一、高斯消元的原理

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

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 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
",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