高斯消元

高斯消元是线性代数规划中的一个算法，可用来为线性方程组求解，高斯消元法可以用在电脑中来解决数千条等式及未知数。

ps：若要解出 \(n\) 个未知数的话，则需要 \(n\) 个有意义的方程。

例如有 \(n\) 个方程组，其中一个是 \(0 \times x = 0 \times y\) 你会发现无论 \(x\) 和 \(y\) 取何值方程都相等，这种方程则无解。

矩阵

高斯消元算法需要用到一个叫矩阵的东西。

有一个 \(A\) 矩阵是系数矩阵，则是记录每一个未知数的系数。

例如：

同时还有一个答案矩阵和一个未知数矩阵。

矩阵乘法

两个大小分别为 \(m \times n\) 和 \(n \times p\) 的矩阵 \(A, B\) 相乘的结果为一个大小为 \(m \times p\) 的矩阵。将结果矩阵记作 \(C\)，则

特殊地，定义 \(A^0\) 为单位矩阵 \(I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}\)。

ps:单位矩阵乘任何矩阵都等于它本身。

算法

好了说了这么多也终于要开始讲高斯消元了。

在高斯消元的矩阵中，因为有 \(n\) 个未知数所以这个矩阵的系数部分是 \(n \times n\) 的。

特殊的在这里我们把答案矩阵与系数矩阵合并在一起，把答案矩阵放在第 \(n + 1\) 列。

了解完初始矩阵后我们开始消元，简单来说消元的过程就是把原矩阵转化为单位矩阵的过程。

----------------------------------------------------------------------------------我是分割线----------------------------------------------------------------------

step1

我们先从当前第 \(i\) 行往第 \(n\) 行找 \(a_i\) 最大的一行并把它替换到当前这行。

for(int j=i;j<=n;j++)
{
	if(abs(a[j][i])>maxx)
	{
		maxx=abs(a[j][i]);
		l=j;
	}
} 
for(int k=1;k<=n+1;k++)
{
	swap(a[i][k],a[l][k]);
}
if(a[i][i]==0)//在这里我们需要判断误解的情况，即最大为0
{
	cout<<"No Solution"<<endl;
	return 0;
}

step2

首先我们需要把第 \(i\) 行的第 \(i\) 个的系数变为一，所以我们需要把这一行所有都除以一个 \(a_i\) 。

for(int j=i+1;j<=n+1;++j)
{
	a[i][j]/=a[i][i];
} 

step3

下一步我们需要把第 \(i + 1\) 行到第 \(n\)行的第 \(a_i\) 个系数都变为0。

所以我们需要把第 \(j\) 行减一个第 \(i\) 行乘第 \(i\) 列。

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

step4

到这里已经快完成了，我们只需要将最后一个答案一次带入求解就行了。

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

code

#include<bits/stdc++.h>
using namespace std;
double max(double a,double b)
{
	if(a>=b) return a;
	if(a<b) return b;
}
int n;
double a[1010][1010];
double a1[1010][1010];
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 maxx=-INT_MAX,l;
		for(int j=i;j<=n;j++)
		{
			if(abs(a[j][i])>maxx)
			{
				maxx=abs(a[j][i]);
				l=j;
			}
		} 
		for(int k=1;k<=n+1;k++)
		{
			swap(a[i][k],a[l][k]);
		}
		if(a[i][i]==0)
		{
			cout<<"No Solution"<<endl;
			return 0;
		}
		for(int j=i+1;j<=n+1;++j)
		{
			a[i][j]/=a[i][i];
		} 
		a[i][i]=1;
		for(int j=i+1;j<=n;j++)
		{
			a1[j][i]=a[j][i];
		} 
		for(int j=i+1;j<=n;++j)
		{
            for(int k=i;k<=n+1;++k)
			{
            	a[j][k]-=a[i][k]*a1[j][i];
			} 
        }
	}
	for(int i=n-1;i>=1;i--)
	{
		for(int j=i+1;j<=n;j++)
		{
			a[i][n+1]-=a[i][j]*a[j][n+1];
		}
	}
    for(int i=1;i<=n;++i)
	{
		printf("%.2f\n",a[i][n+1]);
	} 
	return 0;
}