一、高斯消元的原理
对于n元的m个线性方程组成的方程组,我们将其以矩阵的形式记录下来:
a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn
然后进行初等行列变换,尝试构造出一个上三角矩阵,逐步使系数不为零的项减少;
等最后只剩下一个系数不为零时,进行回代,逐步求出已知解。(详解过程咨询小学老师)
二、高斯消元的实现
老老实实的回代代码参见其他人的博客,这里介绍一种比较毒瘤的不回代暴力消元法:
1.Process
对于每个方程,按照一定的规则(后话)挑选一个主元,记录该主元对应第几个方程,然后用初等行列变换消去其他所有该元的系数;
最后我们得到的是一个每行只有一个数不为零,每列只有一个数不为零的鬼畜矩阵(自己脑补)
此时令ans向量对应的数字出去该行的非零系数,即为对应该元的解。
2.Code
设c数组为已经记录系数的数组(格式见上方),那么c应该是n行n+1列的,最后一列代表方程的常数项;
设w数组记录每个方程的主元是第几项,v数组记录答案;
当多解时输出“Multiple solutions”,无解时输出”No solution”;
double c[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(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(i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][where[i]];
}
3.notice
(1)精度的设置
众所周知浮点数是有精度丢失的,在高斯消元中,精度的选择要依题目而定,精度过低会导致系数较小的数被误认为系数为零,而精度过高也有可能会导致误差大于精度,导致本应该系数消为0的项误认为系数不为零,所以精度的选择是很哲学的问题,要注意。
推荐范围:1e-4到1e-10
(2)主元的选取原则
接着(1)说,我们知道,用浮点数是有精度丢失的,那么让一个较大的数除以一个极小的数误差自然大的可想而知,所以我们想得到在精度允许的条件下系数最大的主元,所以对于每个方程,我们都应该选择最大系数的元作为主元。
(3)在模2意义下的高斯消元
详见相关应用中第三个例题的讲解;