->知识点整理-高斯消元<-
典型题目:
XJOI 1822:Civilization
知识点:
高斯消元其实在小学初中解多元一次方程的时候已经接触过了。其实,高斯消元就是建立在方程中加减消元和乘除消元之上的。只不过,高斯消元法把这两种方法应用于矩阵之中,使得高斯消元的复杂度达到O(n³)(相比于真正的去解方程可是要快的多了,想一想你手解100000元一次方程需要多久就行了)。
然后就是高斯消元的实现了。首先,我们要把高斯消元的每一个系数都加入到N*(N+1)的矩阵中(第N行第N+1列填的是这个方程的答案,但在最后要被我们转化为第N元的解)。
举个例子:假设我们要解如下一个三元一次方程组。
首先,把他们的系数都加入一个矩阵当中:
然后,我们要确定我们所期望的解决后的矩阵是怎样的。这样的矩阵无非两种,一种是上三角矩阵,第i行的元素个数为n-i+1个,并且第i行的前i-1个元素都必须为0。然后从最后一行开始解方程,不断地代入上一行的方程,知道解除答案。而另一种方法则是第i行只有第i个元素不为0,这样可以直接解出每个数的解。然而前者实现较为简单,所以我们今天就用前者。
然后开始我们的消元。首先先从第一列开始,先把该列中绝对值最大的元素所在的哪一行提至第一列(目的是为了提高数据的稳定性,数学运算都是含舍入误差的计算,选主元进行消去可以极大降低舍入误差)。然后就想简单的乘除加减消元一样,把下面的每一列的值都消为0。
下面是第一步举例:
最后是结果的矩阵:
现在只需要从下往上一个一个的解出每一元的解即可。最后得出的解为x=2,y=5,z=8。
然后该怎么判断方程是否有解呢?
只需要按下面的方法判断即可:
① 无解
当方程中出现(0, 0, …, 0, a)的形式,且a != 0时,说明是无解的。
② 唯一解
条件是k = equ,即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。
③ 无穷解。
条件是k < equ,即不能形成严格的上三角形,自由变元的个数即为equ – k。
代码实现:
#include<bits/stdc++.h>
using namespace std;
typedef long double ld;//这里用了long double,是为了精度,但是XJOI原题需要用另一种方法来记录消元时的数,有兴趣可自己去试试
const int maxn=205;
ld a[maxn][maxn];
char in[maxn];//由于XJOI好像不能读long double,所以用了读字符串的方式
void guass(ld x[maxn][maxn],int n)
{
int r;
for(int i=0;i<n;i++)
{
r=i;
for(int j=i+1;j<n;j++)
{
if(fabs(x[j][i])>fabs(x[r][i]))r=j;
}
if(r!=i)for(int j=0;j<=n;j++)swap(x[r][j],x[i][j]);//找出绝对值最大的那一列并进行交换
for(int j=n;j>=i;j--)
{
for(int k=i+1;k<n;k++)
{
x[k][j]-=x[k][i]/x[i][i]*x[i][j];//如果中间用一个中间值进行计算,有可能会造成精度误差,所以我直接进行了计算
}
}
}
for(int i=n-1;i>=0;i--)
{
for(int j=i+1;j<n;j++)
{
x[i][n]-=x[j][n]*x[i][j];
}
x[i][n]/=x[i][i];
}
}
int main()
{
int n;
scanf("%d",&n);
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
{
scanf("%s",in);
int r=0;
for(int k=strlen(in)-1;k>=0;k--)
{
a[i][j]+=(ld)(in[k]-'0')*pow(10,r);
r++;
}
}
}
guass(a,n);
double res;
for(int i=0;i<n;i++)
{
res=(double)a[i][n];//硬生生转回了double输出
printf("%.0lf\n",res);
}
}