高斯消元

高斯消元模板
我们可以这样理解高斯消元：
如：解一个三元方程:

x1+2*x2-x3=6
2*x1+x2-3*x3=9
-x1-x2+2*x3=7

我们可以将其转化成矩阵的形式：

1 2 -1				x1				-6
2 1 -3       *      x2		=  		-9	
-1 -1 2				x3				 7

既然x1，x2，x3是我们要求的未知数，不妨将它们先行舍去，同时将等式右边的数加入到左侧的矩阵内，这就是增广矩阵：

1 2 -1 -6	I
2 1 -3 -9	II
-1 -1 2 7	III

接下来，要解出这几个方程，我们先要了解矩阵的一些基本操作：
初等行（列）变换：
1.交换两行
2.把某一行i的所有元素都乘以k
3.把某一行i的所有元素对应加到j行的所有元素
接下来，我们可以明确我们的目标：
既然该矩阵表示的是

x1+2*x2-x3=6
2*x1+x2-3*x3=9
-x1-x2+2*x3=7

所以当我们把该矩阵化简成如下矩阵(k1，k2，k3等为常数)：

k1	k2	k3	k4
0	k5	k6	k7
0	0	k8	k9

易知该矩阵表示如下等式：

k1*x1+k2*x2+k3*x3=k4
0*x1+k4*x2+k5*x3=k7
0*x1+0*x2+k8*x3=k9

从第三个式子中我们可解出x3
再推到第二个式子我们可解出x2
再推到第一个式子我们可解出x1
观察该矩阵中的“0”
可发现这呈阶梯状，它的名字为阶梯矩阵。
所以我们可以逐渐这么推：

1 2 -1 -6
2 1 -3 -9
-1 -1 2 7

1 2 -1 6
0 -3 -1 3
0 1 1 1

1 2 -1 -6
0 -3 -1 3
0 0 2/3 2

其实在这里我们就可以解出x1，x2，x3的具体值了，但我们可以进一步推成更加理想的形式：

1 2 -1 -6
0 1 1/3 -1
0 0 1 3

1 2 -1 -6
0 1 0 -2
0 0 1 3

1 0 0 13
0 1 0 -2
0 0 1 3

不难发现，这已经变成了一个单位矩阵（除去最右边一列）
综上所述，解多元一次方程组的步骤是：
方程组=>增广矩阵=>阶梯矩阵（初等行变换）=>单位矩阵
此外，可以根据Gauss-Jordan的一个优化：
如这组数据:

x1+2*x2-x3=6
2*x1+x2-3*x3=9
-x1-x2+2*x3=7

我们发现，第一个式子中x1的系数小于第二个式子中x1的系数。所以我们可以交换第一行和第二行，来避免整数除以小数中的精度误差
根据这些步骤，我们可以写出如下代码:

#include <bits/stdc++.h>
using namespace std;
const int maxn=100+10;
double a[maxn][maxn];
int n;
int main(){
	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 maxr=i;
		for(int j=i+1;j<=n;j++)
			if(abs(a[j][i])>abs(a[maxr][i]))
				maxr=j;
		for(int j=1;j<=n+1;j++)
			swap(a[i][j],a[maxr][j]);
		if(a[i][i]==0){
			printf("No Solution\n");
			return 0;
		}
		for(int j=1;j<=n;j++)
			if(j!=i){
				double rate=a[j][i]*1.0/a[i][i];
				for(int k=i+1;k<=n+1;k++)
					a[j][k]-=a[i][k]*rate;
			}
	}
	for(int i=1;i<=n;i++)
		printf("%.2lf\n",a[i][n+1]*1.0/a[i][i]);
	return 0;
}

训练题:
代码：

#include <bits/stdc++.h>
using namespace std;
const int maxn=1000+10;
double a[maxn][maxn],f[maxn][maxn];
int n;
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n+1;i++)
		for(int j=1;j<=n;j++)	
			scanf("%lf",&f[i][j]);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++){
			a[i][j]=2*(f[i][j]-f[i+1][j]);
			a[i][n+1]+=f[i][j]*f[i][j]-f[i+1][j]*f[i+1][j];
		}
	for(int i=1;i<=n;i++){
		int maxr=i;
		for(int j=i+1;j<=n;j++)
			if(abs(a[j][i])>abs(a[maxr][i]))
				maxr=j;
		for(int j=1;j<=n+1;j++)
			swap(a[i][j],a[maxr][j]);
		if(a[i][i]==0){
			printf("No Solution\n");
			return 0;
		}
		for(int j=1;j<=n;j++)
			if(j!=i){
				double rate=a[j][i]*1.0/a[i][i];
				for(int k=i+1;k<=n+1;k++)
					a[j][k]-=a[i][k]*rate;
			}
	}
	for(int i=1;i<=n;i++)
		printf("%.3lf ",a[i][n+1]*1.0/a[i][i]);
	return 0;
}