P3389 【模板】高斯消元法
高斯消元
高斯消元是线性代数规划中的一个算法，可用来为线性方程组求解，高斯消元法可以用在电脑中来解决数千条等式及未知数。
ps：若要解出 \(n\) 个未知数的话，则需要 \(n\) 个有意义的方程。
例如有 \(n\) 个方程组，其中一个是 \(0 \times x = 0 \times y\) 你会发现无论 \(x\) 和 \(y\) 取何值方程都相等，这种方程则无解。
矩阵
高斯消元算法需要用到一个叫矩阵的东西。
有一个 \(A\) 矩阵是系数矩阵，则是记录每一个未知数的系数。
例如：
\[A = \begin{bmatrix} a_{1 1} & a_{1 2} & \cdots & a_{1 n} \\ a_{2 1} & a_{2 2} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m 1} & a_{m 2} & \cdots & a_{m n} \end{bmatrix} \text{.}
\]
同时还有一个答案矩阵和一个未知数矩阵。
矩阵乘法
两个大小分别为 \(m \times n\) 和 \(n \times p\) 的矩阵 \(A, B\) 相乘的结果为一个大小为 \(m \times p\) 的矩阵。将结果矩阵记作 \(C\)，则
\[c_{i j} = \sum_{k = 1}^{n} a_{i k} b_{k j}
\]
特殊地，定义 \(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;
}