高斯消元法——operator_

1.引入

呃，二元一次方程组解过吧，学校教的是加减消元法和代入消元法，代入消元法是要对原方程移项的，有点麻烦，而加减消元法虽然有很多巧妙的解法，但电脑不知道，所以用最一般的方法：加减消元法正推，代入消元法逆向回带。这就是高斯消元法。

2.高斯消元

好像没什么好说的$(?)$纯纯模拟，直接放代码了。

int n;
double a[105][105],x[105];//a[i][j]指的是第i个方程中第j个元素的系数
bool Gauss() {
	for(int i=1;i<=n;i++) {
		int j=i;
		while(j<=n&&fabs(a[j][i])<=eps) j++;//找一个该未知数系数不为0的方程
		if(j>n) return false;//............1
		swap(a[i],a[j]);
		for(j=i+1;j<=n;j++) {//对接下来的每个方程分别消元
			if(fabs(a[j][i])>eps) {
				double tmp=a[j][i]/a[i][i];
				for(int k=i;k<=n+1;k++)
					a[j][k]-=a[i][k]*tmp;
			}
		}
	}
	for(int i=n;i>=1;i--) {
		x[i]=a[i][n+1];
		for(int j=i+1;j<=n;j++)
			x[i]-=x[j]*a[i][j];
		x[i]/=a[i][i];
	}
	return true;
}

下面是返回无解或无数解情况的模板

int Gauss() {
	int r=1;
	for(int i=1;i<=n;i++) {
		int j=r;
		while(j<=n&&fabs(a[j][i])<=eps) j++;
		if(j>n) continue;//............2
		swap(a[r],a[j]);
		for(j=r+1;j<=n;j++) {
			if(fabs(a[j][i])>eps) {
				double tmp=a[j][i]/a[r][i];
				for(int k=i;k<=n+1;k++)
					a[j][k]-=a[r][k]*tmp;
			}
		}
		r++;
	}
	if(r<=n) {
		int flag=0;
		for(int i=r;i<=n;i++) 
			if(fabs(a[i][n+1])>eps) flag=-1;
		return flag; 
	}
	for(int i=n;i>=1;i--) {
		x[i]=a[i][n+1];
		for(int j=i+1;j<=n;j++)
			x[i]-=x[j]*a[i][j];
		x[i]/=a[i][i];
	}
	return 1;
}

差别就在第一个版本标$1$处遇到无唯一解直接返回$false$，第二个版本先跳过，最后再处理是错因

3.高斯-约旦消元法

事实上上面的$2$个版本还有一个优化版本，这就是高斯-约旦消元法，它有两个优化。

优化$1:$高斯消元法的精度误差很大，所以尝试降低误差。你可以尝试一下，会发现除法的误差最大，同时除数越大误差越小 （感性理解）。所以我们在$1$处不再是找到一个就停止，而是找系数最大的方程。

优化$2:$这是使代码更好写的优化。朴素高斯消元法的回带过程能不能优化？当然可以。之所以要回带，是因为高斯消元后的方程组呈三角形，而利用高斯-约旦消元法，我们可以直接得到$n$个形如$kx=b$的简单方程。

代码如下

bool Gauss() {
	for(int i=1;i<=n;i++) {
		int j=i;
		for(int k=i+1;k<=n;k++)
			if(fabs(a[j][i])<fabs(a[k][i])) j=k;
		if(fabs(a[j][i])<eps) return false;
		swap(a[i],a[j]);
		for(j=1;j<=n;j++) {
			if(i!=j) {
				double tmp=a[j][i]/a[i][i];
				for(int k=i;k<=n+1;k++)
					a[j][k]-=a[i][k]*tmp;
			}
        }
	}
	for(int i=1;i<=n;i++)
		x[i]=a[i][n+1]/a[i][i];
	return true;
}