高斯消元专题
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了高斯消元专题相关的知识，希望对你有一定的参考价值。
[高斯消元专题]
最近学了高斯消元，理解起来非常滴简单，就是像我们平常解方程组一样的加减消元，代入消元什么的。但是为了使得这个算法变的比较好梳理，能按照一定的套路去打，那就简单多了。
那么，为了方便得出算法的写法与最终求解的写法，我们一般会将一个方程组一一对应转为矩阵（转为矩阵是什么意思？就是把每一个方程的每一项的系数取下来（包括最后的常数项），然后一一对应的放置在一个矩阵里面）后变成下述形式：
显然，如果某些方程组（矩阵）呈现出这样的形式的话，我们就非常容易能把它解出来。
但是如果初始的方程组不是这样的形式，那怎么办？
还记得加减消元和代入消元不？这样两个方程之间至少能消去一个未知数。
那么，我们做n^3次，连续的这样做消元（当然也要有顺序和策略，这是消元的核心），就能形成上述形式的矩阵。
当然，我们一般选择消成上三角形的矩阵，效率更高更好写。
形成这样的形式后，我们然后从下向上一层一层回代就能求出所以未知数的值了（当然还有无解和多解的情况）。
那么，如何得到一个上三角矩阵呢？每次仅对当前行下面的元素消元。如样例（用最小公倍数法）：
[2] 3 1 16
[1] -2 -2 -7
0 -1 7 18
| |
| |
\\/ \\/
（原1式-2*原2式=新2式）
2 3 1 16
0 [7] 3 30
0 [-1] 7 18
| |
| |
\\/ \\/
（原2式+7*原3式=新3式）
2 3 1 16
0 7 3 30
0 0 52 156
得出x(3)=3，再遍历第二行，将x(3)=3代入，得7x(2)+3*3=30，解出x(2)=3，再代入第一行，得2*x(1)+3*3+3=16，x(1)=2。
下面是几个与方程有关的问题（转自“诚叙”的cnblogs）
Q1：方程无解怎么判断？
A1：如果有一行方程的所有变量的系数都为0，而常数项不为0，方程当然就无解啦。
Q2：方程多解是什么状况呢？
A2：如果有一行方程的所有变量的系数都为0，常数项也为0，那么这就说明出现了自由元（就是在这个方程中不受约束可以自由取值的未知数），且自由元的个数=全部为0的行数，这样就会导致方程多解了。
Q3：方程唯一解是什么样子呢？
A3：不是上面两个不就是了么.....好吧，其实表现在矩阵上是完整美好标准的上三角矩阵或标准型矩阵。
Q4：您上面都说的是解线性的实数或正整数解的方程，但是在题目中经常碰到解异或方程组，这该怎么解呢？
A4：是啊是啊~异或方程组是经常可以在题目中见到的，例如很经典的开关问题....使用的思想还是我们上面讨论的加减消元的思想，转化的矩阵也是我们上面的矩阵，只是矩阵中的系数通常只有0和1了，表示的是如果改变列元素是否会对行元素有影响，例如打开开关1可以让灯2和灯3改变状态，那么a(1,2)=a(1,3)=1，这样就可以通过2，3的状态，求出是不是改变了开关1。（可能还是不太懂....就是我们列的方程成了已知灯的前后状态，求解这些开关是怎么用的，那么我们的已知量就是灯与开关的关系（系数）&灯的状态（常数项），未知量是开关是否使用，大致是一个a(1,1)*x[1] xor a(1,2)*x[2] xor...xor a(1,n)*x[n] 来求解x[1..n]的过程，还是看不懂就看我的另一篇做题的博客好啦）
上面大概介绍了下方程的建立方法，那么我们在解方程的时候我们的消元方法就要变成异或消元了，利用的是1 xor 1=0，如果你要解某个变量 i 就找到一行 i 的系数不为0的交换到当前行（如果当前行都没有这个元素，就解不了了（上面普通方程求解的时候也可以这样）），然后找到其他当前变量的系数不为0的，异或之（每个系数都异或，包括常数），最后也能化成标准或上三角矩阵，从而求解。
还要特殊说明的是在异或方程中如果知道了方程的自由元数量就可以知道总共有多少组解啦！因为只有1或0两种可能嘛，所以答案就是2^n(n为自由元的数量)
Q5：高斯消元还有什么用呢？
A5：还可以和期望概率扯上关系呢，因为一个事件的期望值很可能和与它相关的其他事件扯上关系，例如我从原点出发，每次只能走一步，往上或往右，那么走到(i,j)的期望步数就是从(i-1,j)走来的期望步数*从(i-1,j)走来的概率+1再加上从(i,j-1)走来的期望步数*从(i,j-1)走来的概率+1，是不是很像方程组呢.....通过这样也能建立方程组（....啊啊啊其实我开始有点口糊的感觉了....我还没打过这方面的题~加油你们去膜大神们吧~moto_No.1还讲了个在线求高斯，大家可以去看看他的博客，概率什么的，查查题目就能找到啦(有趣的“驱赶猪猡”....)）
这位dalao写的很不错哦~~~
那么，下面就给出几份不同情况的代码：
对于普通的整数方程组（非实数）：
1 int Gauss(int equ,int var){ 2 int R=1,C=1; 3 for (; R<=equ&&C<=var; R++,C++){ 4 int Mx=R; 5 for (int i=R+1; i<=equ; i++) 6 if (abso(a[i][C])>abso(a[Mx][C])) i=Mx; 7 if (Mx!=R) 8 for (int j=C; j<=var+1; j++) swap(a[Mx][j],a[R][j]); 9 if (a[R][C]==0){R--; continue;} 10 for (int i=R+1; i<=equ; i++) if (a[i][C]!=0){ 11 int LCM=lcm(abso(a[i][C]),abso(a[R][C])); 12 int tx=LCM/abso(a[i][C]),ty=LCM/abso(a[R][C]); 13 if (a[i][C]*a[R][C]<0) ty=-ty; 14 for (int j=C; j<=var+1; j++) a[i][j]=a[i][j]*tx-a[R][j]*ty; 15 } 16 } 17 for (int i=R; i<=equ; i++) if (a[i][var+1]!=0) return -1; //无解的情况 18 if (R<=var) return var-R+1; //多解的情况 19 for (int i=n; i; i--){ 20 int res=a[i][var+1]; 21 for (int j=i+1; j<=var; j++) if (a[i][j]!=0) res-=a[i][j]*x[j]; 22 x[i]=res/a[i][i]; 23 } 24 return 0; //唯一解的情况 25 }
对于某些01异或方程组（给出的是有唯一解的情况，多解和无解想必也不难）：
1 void Gauss(){ 2 int row=1,col=1; 3 for (; row<=N&&col<=N; row++,col++){ 4 int Mx_r=row; 5 for (int i=row+1; i<=N; i++) 6 if (abso(a[i][col])>abso(a[Mx_r][col])) Mx_r=i; 7 if (Mx_r!=row) 8 for (int j=col; j<=N+1; j++) swap(a[Mx_r][j],a[row][j]); 9 for (int i=row+1; i<=N; i++) if (a[i][col]) 10 for (int j=col; j<=N+1; j++) a[i][j]^=a[row][j]; 11 } 12 for (int i=N; i>=1; i--){ 13 x[i]=a[i][N+1]; 14 for (int j=i+1; j<=N; j++) x[i]^=(a[i][j]&&x[j]); 15 } 16 }
对于某些同余方程组（模数必须为质数）：
1 int Gauss(int equ,int var){ 2 int R=1,C=1; 3 for (; R<=equ&&C<=var; R++,C++){ 4 int Mx=R; 5 for (int i=R+1; i<=equ; i++) 6 if (abso(a[i][C])>abso(a[Mx][C])) Mx=i; 7 if (Mx!=R) 8 for (int j=C; j<=var+1; j++) swap(a[Mx][j],a[R][j]); 9 if (a[R][C]==0){R--; continue;} 10 for (int i=R+1; i<=equ; i++) if (a[i][C]!=0){ 11 int LCM=lcm(abso(a[i][C]),abso(a[R][C])); 12 int tx=LCM/abso(a[i][C]),ty=LCM/abso(a[R][C]); 13 if (a[i][C]*a[R][C]<0) ty=-ty; 14 for (int j=C; j<=var+1; j++) 15 a[i][j]=(((a[i][j]*tx-a[R][j]*ty)%TT)+TT)%TT; 16 } 17 } 18 for (int i=R; i<=equ; i++) if (a[i][var+1]!=0) return -1; 19 if (R<=var) return var-R+1; 20 for (int i=n; i; i--){ 21 int res=a[i][var+1]; 22 for (int j=i+1; j<=var; j++) 23 if (a[i][j]!=0) res-=a[i][j]*x[j],res=((res%TT)+TT)%TT; 24 x[i]=((res*Q_pow(a[i][i],TT-2))%TT+TT)%TT;// x[i]=res/a[i][i]求逆元 25 } 26 return 0; 27 }
对于某些实数方程组（唯一解情况）：
1 void Gauss(int equ,int var){ 2 int R=1,C=1; 3 for (; R<=equ&&C<=var; R++,C++){ 4 int Mx=R; 5 for (int i=R+1; i<=equ; i++) 6 if (dcmp(dabs(a[i][C])-dabs(a[Mx][C]))==1) Mx=i; 7 if (Mx!=R) 8 for (int j=C; j<=var+1; j++) swap(a[Mx][j],a[R][j]); 9 if (dcmp(a[R][C])==0){R--; continue;} 10 for (int i=R+1; i<=equ; i++) if (dcmp(a[i][C])!=0){ 11 double k=a[i][C]/a[R][C]; 12 for (int j=C; j<=var+1; j++) a[i][j]=a[i][j]-k*a[R][j]; 13 } 14 } 15 for (int i=var; i; i--){ 16 double res=a[i][var+1]; 17 for (int j=i+1; j<=var; j++) res-=a[i][j]*x[j]; 18 x[i]=res/a[i][i]; 19 } 20 }
对于某些实数方程组多解的情况（求出哪些是自由元，以及非自由元的值）：
1 int Gauss(int equ,int var){ 2 int R=1,C=1; 3 for (; R<=equ&&C<=var; R++,C++){ 4 int Mx=R; 5 for (int i=R+1; i<=equ; i++) 6 if (dcmp(dabs(a[i][C])-dabs(a[Mx][C]))==1) Mx=i; 7 if (Mx!=R) 8 for (int j=C; j<=var+1; j++) swap(a[Mx][j],a[R][j]); 9 if (dcmp(a[R][C])==0){R--; continue;} 10 for (int i=R+1; i<=equ; i++) if (dcmp(a[i][C])!=0){ 11 double k=a[i][C]/a[R][C]; 12 for (int j=C; j<=var+1; j++) 13 a[i][j]=a[i][j]-k*a[R][j]; 14 } 15 } 16 for (int i=R; i<=equ; i++) if (dcmp(a[i][var+1])!=0) return -1; 17 if (R<=var){ 18 for (int i=R-1; i; i--){ 19 int num=0,idx; 20 for (int j=1; j<=var; j++) if (freee[j]&&dcmp(a[i][j])!=0){ 21 num++,idx=j; 22 } 23 if (num!=1) continue; 24 double res=a[i][var+1]; 25 for (int j=1; j<=var; j++) if (j!=idx) 26 if (dcmp(a[i][j])!=0) res-=a[i][j]*x[j]; 27 x[idx]=res/a[i][idx],freee[idx]=0; 28 } 29 return var-R+1; 30 } 31 for (int i=1; i<=var; i++) freee[i]=0; 32 for (int i=var; i; i--){ 33 double res=a[i][var+1]; 34 for (int j=i+1; j<=var; j++) res-=a[i][j]*x[j]; 35 x[i]=res/a[i][i]; 36 } 37 return 0; 38 }
当然还有其他的情况，那代码其实并不重要了，重要的是我们要能够举一反三。->_->
下面是习题：
以上是关于高斯消元专题的主要内容，如果未能解决你的问题，请参考以下文章