ATP记得它在很久以前看过一点点高斯消元的东西然后做过一点点题目。。但是当时实在是太zz了所以本来就没有很懂这个东西现在更是忘得差不多了。。
所以现在就当重新学一遍了QwQ
一点口胡的解释
高斯消元。。听起来这个名字很高大上实际上它也确实很高大上但是它的操作过程实际上和数学课学的加减消元法是一个道理的啦。
它可以用来求解线性方程组 Ax=b ,其中 A 是一个
高斯消元还可以求解可逆矩阵的逆矩阵。还可以求解矩阵的秩。(秩是啥?我不会= =)
高斯消元的操作过程实际上非常暴力。。总的来说分两个步骤,首先要通过矩阵行变换也就是加减消元消来消去把这些方程变成某个等价形式。这个等价形式的特征体现在系数矩阵A里就是A变成了一个“上三角矩阵”,也就是只有左上-右下对角线及其上方有数字,这条对角线的下面都是0,大概就是这个样的东西:
这样的话可以看出来第i个方程有n-i+1个系数,那么最后一个方程只有一个系数,也就是经过这样消元以后最后一个方程变成了 An,n∗xn=bn 的形式,那么就可以直接得到 xn 的值了。得到 xn 以后就可以往回代,因为倒数第2个方程有2个系数,其中一个系数是关于 xn 的,把得到的 xn 的值带回去以后就可以把它消成只有一个未知量的方程,那么就可以得到 xn−1 。同理把已经求出的未知量都往回代,就可以求出所有的未知数。
贴一个板子。出自BZOJ1013球形空间产生器,是高斯消元的一个板子题:
void Gauss_Eli(){
int num;
for (int i=1;i<n;i++){
num=i;
for (int j=i+1;j<=n;j++)
if (fabs(A[j][i])>fabs(A[num][i]))
num=j;
for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);
swap(b[num],b[i]);
for (int j=i+1;j<=n;j++){
if (fabs(A[j][i])<=eps) continue;
double t=A[j][i]/A[i][i];
for (int k=1;k<=n;k++)
A[j][k]-=A[i][k]*t;
b[j]-=b[i]*t;
}
}
for (int i=n;i>=1;i--){
ans[i]=b[i]/A[i][i];
for (int j=i;j>=1;j--)
b[j]-=A[j][i]*ans[i];
}
}
可以看出高斯消元具体的实现复杂度是 O(n3) 的。这里以及下面我们都把将要留在上三角矩阵的第i行的方程称为未知数i的关键方程,因为如果未知数 xi 的值是唯一的,上三角矩阵的第i行必须有未知数 xi 的系数。
上面的代码大体来说就是先枚举每一个变量i找到它的关键方程,然后把它的关键方程下面所有方程 xi 的系数都消成0让它符合上三角矩阵的形式,最后再进行回代。
首先我们要从上到下枚举每一行。设当前枚举到的是第i行,那么前i-1行已经消成了上三角矩阵的一部分,不能再动了。我们要找第i个未知量的关键方程只能从后面找。上面代码的第4..9行干了一个看起来比较奇怪的事情:从后面的方程中找了一个未知数i的系数最大的方程然后换到第i个位置。据说找一个最大的能减少浮点误差?然而比较实质性的作用还是保证选中的这个“准”关键方程里 xi 的系数不是0吧。。如果它题目里给的方程排列比较凑巧,第i行上 xi 的系数正好是0,那么它肯定不能作为 xi 的关键方程,肯定要找一个系数不是零的搞过来用。所以这一步还是比较重要的。
找到了“准”关键方程以后我们要把它变成关键方程。要干的事情就是让它下面的所有方程 xi 的系数都变成0。那么就要用到加减消元法了。设第i行下面的某一行为j,如果第j行中 xi 的系数不为0,那么我们就要把它变成0。联想一下数学课上解二元一次方程组的方法,我们要把某一个方程中某一个未知量消掉,方法就是把这个方程扩大一个特定的倍数,让要消掉的这个未知量的系数和另一个方程的对应系数相等,再把这两个方程相减就能消掉目标未知数。那么这里目标未知数 xi 的系数是 Ai,i ,我们要把它扩大某个倍数让 Ai,i 变成 A