浅谈高斯消元

引入

先来看一个一元一次方程：
$$3x-9=0$$
对于这样一个形如$ax+b=0$的方程，我们很容易的解得$x=3$。

在电脑程序中，我们的处理也很简单：

double solve(double a,double b){  //ax=b
	return b/a;
}

我们再看一个二元一次方程组：
$$\{\begin{array}{r}2x+5y=9①\\ 4x-3y=5②\end{array}$$
对于人脑来讲，我们可以一眼看出该方程的解为$x=2,y=1$，但是电脑不行。

为了让这种解法适用于电脑，我们需要还原我们的思考过程。

数学书上教我们的方法是加减消元法：

$①\times 3+②\times 5得,26x=52$

$解得x=2$

$将x=2代入②式得8-3y=5$

$解得y=1$

由上面的过程我们可以看出，加减消元法将一个含有多个未知数的方程组转化成了多个简单的一元一次方程。

将上面的做法扩展到$n$元$n$次方程组，就是我们今天要介绍的高斯消元。

高斯消元

原理

高斯消元的原理基于$3$条定理：

$1.$两方程互换，解不变。

$2.$一方程乘上非零数$k$，解不变。

$3.$一方程加上另一方程，解不变。

实现

我们以这道题为例，

我们的输入是：

1 2 1 7
2 -1 3 7
3 1 2 18

变成常见的方程组：
$$\{\begin{array}{r}x+2y+z=7①\\ 2x-y+3z=7②\\ 3x+y+2z=18③\end{array}$$
我们将其转化成一个$n\times (n+1)$的矩阵：
$$\left[\begin{array}{ccccc}1& 2& 1& |& 7\\ 2& -1& 3& |& 7\\ 3& 1& 2& |& 18\end{array}\right]$$

我们再进行除法时会产生精度误差，故我们使$x$最大的在第一行，矩阵变为：

$$\left[\begin{array}{ccccc}3& 1& 2& |& 18\\ 2& -1& 3& |& 7\\ 1& 2& 1& |& 7\end{array}\right]$$

我们设矩阵的$3$行从上到下分别为$R_1,R_2,R_3$，就可以按照常规思路解方程了。

由于方便的考虑，我们需要在消元时把式子里的某一个系数变为$1$。

$R_2\times 3-R_1\times 2=R_2\times 1.5-R_1:$
$$\left[\begin{array}{ccccc}3& 1& 2& |& 18\\ 0& -2.5& 2.5& |& -7.5\\ 1& 2& 1& |& 7\end{array}\right]$$
$R_3\times 3-R_1\times 1:$
$$\left[\begin{array}{ccccc}3& 1& 2& |& 18\\ 0& -2.5& 2.5& |& -7.5\\ 0& 5& 1& |& 3\end{array}\right]$$
$R_3-R_2\times (-2)=R_3\times (-0.5)-R_2:$
$$\left[\begin{array}{ccccc}3& 1& 2& |& 18\\ 0& -2.5& 2.5& |& -7.5\\ 0& 0& -3& |& 6\end{array}\right]$$
转化成方程组：
$$\{\begin{array}{r}3x+y+z=18\\ -2.5y+2.5z=-7.5\\ -3z=6\end{array}$$
然后解出$z$，由$z$解出$y$，再解出$x$，这个操作被称作回代。

这就是加回代的高斯消元解方程的过程。

细节

消除误差

对于浮点数的运算，误差是一个很常见的事情。

我们假设最后解方程的某系数为$0$，但经过多重运算后，当前的系数变为$10^{-10}$，而我们的程序就会判断此数不为$0$，然后运算出一些奇怪的结果。

对此的解决方案就是定义一个极小值$eps$，一般取$10^{-6}$到$10^{-8}$间的数，如果一个数的绝对值小于$eps$,那么我们就认为该数为$0$。

#define eps 1e-6
bool check(double x){
	if(fabs(x)<=eps)  return false;
	return true;
}

判断无解或无穷解

先讨论无解的情况。

从数学角度考虑，只有当$0=非零数$时我们的方程才无解，故如果我们在回代过程中发现消去所有已知量后，某方程的系数全为$0$但常数项非零，那么我们就退出。

int Gauss(){
	/*消元+回代*/
	if(!check(a[i][i])&&check(a[i][n+1]))  return -1;
}

再考虑无穷解的情况。

与上面的情况相似，但当某方程回代后常数项和未知数全为$0$时，方程有无穷解。

int Gauss(){
	/*消元+回代*/
	if(!check(a[i][i])&&!check(a[i][n+1]))  return -1;
}

消元的过程

由于在计算$a{R}_n-b{R}_m$时我们需要将其中一个系数变为$1$，所以我们让方程组上面的式子的系数大于下面的方程组。

int p=i;
for(int j=1;j<=n;j++)
  if(fabs(a[j][i])>fabs(a[p][i]))  p=j;
for(int j=1;j<=n+1;j++)  swap(a[i][j],a[p][j]);

然后我们进行加减消元。

for(int j=i+1;j<=n;j++){
	double tmp=a[j][i]/a[i][i];
	for(int k=1;k<=n+1;k++){
		a[j][k]=a[j][k]-a[i][k]*tmp;
	}
}

回代的过程

将所有已知变量的和减去，判断一下是否有唯一解即可。

for(int i=n;i>=1;i--){
	for(int j=n;j>i;j--){
		if(a[i][j]!=0)  a[i][n+1]-=ans[j]*a[i][j];
	}
	ans[i]=a[i][n+1]/a[i][i];
}

非回代版高斯消元

优点：精度更好，代码更简单（没有回代的过程）

实现

选择一行的某个未知数为主元，直接消掉其余行的所有带该未知数的项，最后得到的答案矩阵就直接为：
$$\left[\begin{array}{ccccccc}{a}_1& 0& 0& ...& 0& |& A_1\\ 0& a_2& 0& ...& 0& |& A_2\\ 0& 0& a_3& ...& 0& |& A_3\\ ...& ...& ...& ...& ...& |& ...\\ 0& 0& 0& ...& a_n& |& A_n\end{array}\right]$$
答案$x_i=A_i/a_i$

for(int j=i+1;j<=n;j++){
	double tmp=a[j][i]/a[i][i];
	for(int k=1;k<=n+1;k++){
		a[j][k]=a[j][k]-a[i][k]*tmp;
	}
}

而我们实现时的具体差异就是将上方代码的

for(int j=i+1;j<=n;j++)

改为

for(int j=1;j<=n;j++)

又由于当我们处理到某一行时，其前面的一些系数已经被前面几行的操作清零了，所以我们枚举$k$的时候只需从$i+1$枚举到$n+1$就可以了。

for(int k=i+1;k<=n+1;k++){
	a[j][k]-=a[i][k]*tmp;
}

Code

#include<bits/stdc++.h>
using namespace std;
double a[105][105];
int main(){
	int n;
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n+1;j++){
			scanf("%lf",&a[i][j]);
		}
	}
	for(int i=1;i<=n;i++){
		int maxn=i;
		for(int j=i+1;j<=n;j++){
			if(fabs(a[j][i])>fabs(a[maxn][i])){
				maxn=j;
			}
		}
		for(int j=1;j<=n+1;j++){
			swap(a[i][j],a[maxn][j]);
		}
		if(!a[i][i]){
			cout<<"No Solution\n";
			return 0;
		}
		for(int j=1;j<=n;j++){
			if(j==i)  continue;
			double tmp=a[j][i]/a[i][i];
			for(int k=i+1;k<=n+1;k++){
				a[j][k]-=a[i][k]*tmp;
			}
		}
	}
	for(int i=1;i<=n;i++){
		printf("%.2lf\n",a[i][n+1]/a[i][i]);
	}
}

注意：所有使用double的地方不要用int！