高斯消元 学习笔记

用于求解方程组。

给定 \(n\) 个关于 \(m\) 个变量的方程组，需要你判断该方程组是否无解、有无数解、有唯一解，并输出唯一的解。

考虑使用消元法。我们枚举一个变量 \(i\) ，从所有没有被操作过的方程式中选出一个，然后用它对其他没有被操作过的方程式进行消元，并将被选中的那个方程式视为被操作过的。每次在消元时我们都要将其余方程式中的变量 \(i\) 的系数变成 0 。处于对精度的考虑，我们会选择系数的绝对值最大的那个方程式进行消元（意会一下：设我们选中的方程式为第 \(k\) 个，我们将该方程式的所有系数全部除 \(i\) 的系数 \(a_{k,i}\) ，而我们将 \(\frac{1}{a_k,i}\) 视为一个基本单位，这个单位越小，我们最后在消元时得到的值损失精度的程度会降低）。

现在出现了一个问题，如果当某一个变量 \(i\) ，所有未操作过的方程式关于它的系数都为 0 怎么办？直接跳过该变量就好了，因为此时剩余的方程式的性质没有变，都是满足前面所有枚举过的变量的系数为 0。同时，这种情况会使得我们在遍历了所有变量后，在矩阵的底部存在一些所有变量的系数均为 0 的方程式，显然，这种方程式的答案一定为 0。如果不为 0，则方程组无解；如果为 0，则方程式中存在一个自由元（可以取任意值），此时方程组有无数解。

将特殊情况排除后，我们看看如何计算一个唯一解。我们发现，对于最后的一个方程式，它的系数只有一个不为 0，于是我们可以算出来。再看上一个方程式，发现其中可能有两个系数不一定为 0，但是其中一个变量已经被我们算出来了，以此类推。于是我们从下往上回带就好了。

Code

#include<bits/stdc++.h>
using namespace std;
const int N=105;const double eps=1e-9;
double a[N][N],ans[N];
int main()
{
	int n,i,j,k,l;double x;
	scanf("%d",&n);
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n+1;j++)
		scanf("%lf",&a[i][j]);
	}
	for(i=1,l=1;i<=n;i++)
	{
		x=abs(a[l][i]);k=l;
		for(j=l+1;j<=n;j++)
		{
			if(abs(a[j][i])>x)
			{
				k=j;x=abs(a[j][i]);
			}
		}
		if(x<=eps)
		continue;
		for(j=i;j<=n+1;j++)
		swap(a[k][j],a[l][j]);
		x=a[l][i];
		for(j=i;j<=n+1;j++)
		a[l][j]/=x;
		for(j=l+1;j<=n;j++)
		{
			for(k=n+1;k>=i;k--)
			a[j][k]=a[j][k]-a[j][i]*a[l][k];
		}
		l++;
	}
	for(i=n;i>=l;i--)
	{
		if(abs(a[i][n+1])>eps)
		{
			printf("-1");return 0;
		}
	}
	if(l<=n)
	{
		printf("0");return 0;
	}
	for(i=n;i>=1;i--)
	{
		for(j=i+1;j<=n;j++)
		a[i][n+1]=a[i][n+1]-a[i][j]*ans[j];
		ans[i]=a[i][n+1]/a[i][i];
	}
	for(i=1;i<=n;i++)
	printf("x%d=%0.2lf\n",i,ans[i]);
	return 0;
}