高斯消元入门

前言

学高斯消元之前，我觉得这东西真难。学完之后，我发现高斯消元其实也挺简单的。

先通过一个例子来初步认识高斯消元

这里有一组方程：

\[\begin{cases}5x+4y=24,\\4x+3y=19.\end{cases} \]

这组方程怎么用高斯消元来求解呢？

我们先用一个\(delta\)来记录要想与②式相加消掉\(x\)，①式所要扩大的倍数。不难发现，在这个例子中，\(delta=-\frac45\)，因此，我们将①式左右两边同乘\(delta\)，得：

\[-4x-\frac{16}5y=-\frac{96}5\tag{等式性质} \]

将②式加上③式，可消掉\(x\)得

\[-\frac15y=-\frac15\tag{加减消元法} \]

结合①式和消元得到的④式，我们就得到了：

\[\begin{cases}5x+4y=24,\\-\frac15y=-\frac15\end{cases} \]

然后，我们就可以从最后一个式子开始，求出每一个未知数的答案：

• ④式：
  \(y=-\frac15÷-\frac15=1\)
  代入①式得

\[5x=20\tag{代入消元法} \]

• ①式：
  \(x=20÷5=4\)

综上所述，原方程组的解为

\[\begin{cases} x=4,\\ y=1. \end{cases} \]

整理思路

接下来，我们整理一下刚才的步骤与思路：

• 第一步：先选择第一行的第一个元素，然后通过加减消元法消去每一行的这一个元素。
• 第二步：我们发现最后一个方程已经可以直接求解，然后用代入消元法将最后一个方程求出来的解代入前面的每一个方程实现消元。

上面的步骤只是对于二元方程的，对于一个\(n\)元方程应该这样用高斯消元法：

1. 新建一个变量\(i\)，初始化\(i=1\)
2. 找到一个行数大于等于\(i\)且第\(i\)个元素系数不为0的方程，将其移至第\(i\)行
3. 如果找不到，就输出\(No\) \(Solution\)，并退出程序
4. 新建变量\(j\)，初始化\(j=i+1\)
5. 用一个变量\(delta\)记录在将想要将第\(i\)行与第\(j\)行相加消掉第\(i\)个元素第\(i\)行所要扩大的倍数，即\(-\frac{第j行第i个元素的系数}{第i行第i个元素的系数}\)
6. 将第\(j\)行每一个元素的系数分别加上第\(i\)行对应的元素的系数乘以\(delta\)后的值
7. 将\(j\)加1，如果\(j\)不大于\(n\)，返回步骤5
8. 将\(i\)加1，如果\(i\)不大于\(n\)，返回步骤2
9. 将\(i\)重新赋值为\(n\)
10. 如果第\(i\)行第\(i\)个元素的系数为0，就输出\(No\) \(Solution\)，并退出程序
11. 否则，该行其他元素的系数肯定都为0，因此第\(i\)个未知数的值为\(\frac{第i行等号右边的值}{第i行第i个元素的系数}\)
12. 将第\(i\)个元素的值代入第\(1\sim i-1\)行的式子中消去第\(i\)个未知数
13. 将\(i\)减1，如果\(i\)不小于1，返回步骤10

再接下来就是代码的实现了。

代码

#include<bits/stdc++.h>
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define abs(x) ((x)<0?-(x):(x))
#define LL long long
#define ull unsigned long long
#define tc() (A==B&&(B=(A=ff)+fread(ff,1,100000,stdin),A==B)?EOF:*A++)
#define N 100
char ff[100000],*A=ff,*B=ff;
using namespace std;
const double eps=1e-6;
int n;
double a[N+5][N+5],s[N+5];
inline int read()
{
	int x=0,f=1;static char ch;
	while(!isdigit(ch=tc())) f=ch^'-'?1:-1;
	while(x=(x<<3)+(x<<1)+ch-48,isdigit(ch=tc()));
	return x*f;
}
inline void swap(double &x,double &y)
{
	double t=x;x=y,y=t;
}
inline void Find_line(int x)//找到一个行数大于等于x且第x个元素系数不为0的方程，将其移至第x行
{
	register int i=x,j;
	while(i<=n&&fabs(a[i][x])<eps) ++i;//只要i不大于n且这一行第x个元素为0，就将i加1
	if(i>n) puts("No Solution"),exit(0);//如果无解，输出"No Solution"，并退出程序
	for(j=1;j<=n;++j) swap(a[x][j],a[i][j]);//将第i行移至第x行
}
int main()
{
	register int i,j,k;
	for(n=read(),i=1;i<=n;s[i++]=read()) for(j=1;j<=n;++j) a[i][j]=read();//读入每一个元素的系数以及等于号右边的值
	for(i=1;i<=n;++i) 
	{
		Find_line(i);
		for(j=i+1;j<=n;++j)//消去[i+1~n]中每一行第i个元素
		{
			double delta=-a[j][i]/a[i][i];//delta记录在将想要将第i行与第j行相加消掉第i个元素第i行所要扩大的倍数
			for(k=i;k<=n;++k) a[j][k]+=a[i][k]*delta;//将第j行每一个元素的系数分别加上第i行对应的元素的系数乘以delta后的值
			s[j]+=s[i]*delta;
		}
	}
	for(i=n;i;--i) 
	{
		if(!a[i][i]) return puts("No Solution"),0;//如果第i行第i个元素的系数为0，就输出"No Solution"，并退出程序
		for(s[i]/=a[i][i],j=i-1;j;--j) s[j]-=a[j][i]*s[i];//计算出第i个未知数的值，并将第i个元素的值代入第1~i-1行的式子中消去第i个未知数
	}
	for(i=1;i<=n;++i) printf("%.2lf\n",s[i]);//输出每一个未知数的值
	return 0;
}

例题

一道例题：【BZOJ1013】[JSOI2008] 球形空间产生器