高斯消元之乱搞矩阵

前言

今日，冰红茶同志跑过来问我高斯消元。他说看了我的博客没有看懂。好吧，这篇博客确实太水了，什么也没讲好，今天重打一发。

—————————————分割线—————————————

2022.3.26 更新：今天又把高斯消元考古了一遍emmm

1.引言

作为一名经常与小学数学打交道的OIER，大家应该都会解多元一次方程组吧~~~

小学老师都讲解过，要想解出包含有多个未知数的方程组，最重要的就是一个个的去消元，在回带。

最后解出方程组的正解。

今天蒟蒻就为大家讲解一下，在遭遇大量的多元一次方程组时的解决方法：

高斯-约旦消元法！！！

要想学会高斯-约旦消元法，就要先了解矩阵；

在列出 ,

$$\left\{ \begin{array}{c} a_1*x1+b_1*x2+c_1*x3=d1\\ a_2*x1+b_2*x2+c_2*x3=d2\\ a_3*x1+b_3*x2+c_3*x3=d3 \end{array} \right.$$

这样的方程组时，我们可以把它转换为一个矩阵：

$$\left\{ \begin{array}{c} a1 , b1 , c2 | d1\\ a2,b2,c2|d2\\ a3,b3,c3|d3 \end{array} \right\}.$$

这样一来，我们就可以利用高斯三板斧：

$$1.用一非零数乘某个线性方程;\\ 2.把一个线性方程的倍数加到另一个线性方程;\\ 3.互换两个线性方程的位置.\\ $$

将其转化为：

$$\left\{ \begin{array}{c} 1,0,0|x1\\ 0,1,0|x2\\ 0,0,1|x3 \end{array} \right\}.$$

就得到了方程的解了~~

2.计算机实现

更：
高斯消元在我们平常的计算中很少运用，但是在计算机求解方程上却非常有用。

因为高斯消元的本质就是三种操作，计算机很好模拟这样的过程。

先讲一下计算机高斯消元的流程：

1：枚举每一列，也就是枚举每一个未知数：

2：先找到这个未知数的主元，我们一般把系数最大的作为主元，因为可以减小误差：

3.将当前的行与主元交换，因为我们必须在当前行进行操作，并求出未知数；

4.在与当前行交换后，从后往前，该未知数后面的所用未知数的系数除以当前主元的系数，这样就对当前行进行了消元。

5.当前行消元后，再从上到下倒着依次消元，注意是倒着消元，而且不是当前行。

6.每一行最后的就是答案了。

3.代码实现

接下来就看代码的实现了：
1.先来一道模板：P3389高斯消元

总体比较简单；
这道题不要求我们考虑无解的情况，所以较为简单。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const int N=1e4;
double A[N][N]; 
int n;
int main()
{
	cin>>n;
	for(int i=0;i<n;i++){
		for(int j=0;j<=n;j++){
			cin>>A[i][j];
		}
	}
	for(int j=0;j<n;j++){//枚举每一列消元
		int i;
		for(i=j;i<n;i++){//找主元
			if(A[i][j]){
				break;
			}
		}
		if(A[i][j]==0){
			puts("No Solution");
			return 0;
		}
		for(int k=0;k<=n;k++){
			swap(A[i][k],A[j][k]);//与当前行交换
		}
		for(int k=n;k>=j;k--){
			A[j][k]/=A[j][j];//消元
		}
		for(int i=0;i<n;i++)
			if(i!=j)
				for(int k=n;k>=j;k--)
					A[i][k]-=A[j][k]*A[i][j];//消去后面每一行的元
	}
	for(int i=0;i<n;i++)
	printf("%.2lf\n",A[i][n]); 
	return 0;
}

2.然后就是我们平常所用的模板了：

点击查看代码

void gauss(){
	for(int i=1;i<=n;i++){
		int max=i;
		for(int j=i+1;j<=n;j++){
			if(fabs(a[j][i])>fabs(a[max][i])) swap(max,j);
		}
		if(max!=i){
			for(int j=1;j<=n+1;j++){
				swap(a[i][j],a[max][j]);
			}
		}
		if(fabs(a[i][i])<eps) continue;
		for(int j=1;j<=n;j++){
			if(j!=i){
				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++){
		a[i][n+1]/=a[i][i];
	}
}

同学们，你们应该都懂了吧~~~

4.例题

[JSOI2008]球形空间产生器

点击查看代码

#include<bits/stdc++.h>
using namespace std;
double a[20][20],A[20][20];
int n;
int main(){
	cin>>n;
	for(int i=0;i<=n;i++){
		for(int j=0;j<n;j++){
			cin>>a[i][j];
		}
	}
	for(int i=0;i<=n;i++){
		for(int j=0;j<n;j++){
			A[i][j]=2*(a[i][j]-a[i+1][j]);
			A[i][n]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
		}
	}
	
	for(int j=0;j<n;j++){
		int i;
		for(i=j;i<=n;i++){
			if(A[i][j]){
				break;
			}
		}
		for(int k=0;k<=n;k++){
			swap(A[i][k],A[j][k]);
		}
		for(int k=n;k>=j;k--){
			A[j][k]/=A[j][j];
		}
		for(int i=0;i<=n;i++){
			if(i!=j){
				for(int k=n;k>=j;k--){
					A[i][k]-=A[j][k]*A[i][j];
				}
			}
		}
	}
	for(int i=0;i<n;i++){
		printf("%.3lf ",A[i][n]);
	}
	return 0;
}