高斯消元法

高斯消元问题通常用来对线性方程组的求解

根据线性代数中所学，线性方程通常可表示成

a[i][0]*x0 + a[i][1]*x1 +a [i][2]*x2 ... + a[i][n-1]*xn-1 = b[i]

0<=i<n

这样就可以得到一个由n个线性方程组成方程组

我们通常将a[][]放入数组中的对应位置，b[] 作为一整列插入到a数组的最后一排

 

 

不断通过对行进行操作

比如 ri - k*rj , 用第i行的元素每一个都减去第j行对应列的元素的k倍，我们总是希望能够通过这个操作不断将矩阵朝着除去 b[]后的上三角

矩阵化简，就是令下三角矩阵均为0

那么就是不断从第一列开始不断将需要化简为0的项转化为0即可

那么我们就可以利用回代的方法，先算出最后一项x，然后不断往前递推，在第i行得到x[i]

 

 

当然这是最好的情况，解决线性方程组往往给的方程的数量和对应x的变量数量不一定相同

我令方程数量为equ ， 变量数量为 var ， 那么如果equ < var 的时候必然是有解的

因为至少有 var-equ+1个自由变量

因为每一行都会解决至多一个自由变量，当然如果某一行和前面的一行 的系数都存在倍数关系 exp: ai = k*aj

那么那一行其实效果是没有的，这个时候就一个自由变量都没有解决了

 

另外如果在化简矩阵的过程中某一行出现了（0,0,0,0....,0,b[i])（b[i] != 0)的情况那么就说明是不存在解的

 

所以对于我们来说，只要除去那些有相同效果的方程，剩下的方程组含有的方程为equ，那么

自由变量的个数就是 var - equ + 1 ， 当自由变量个数大于1的时候我们就可以认为是有无穷多个解的

 

我们利用这个思想解决计算机上的问题时，可以套用如下模版：

 

  1 int gcd(int a , int b)
  2 {
  3     if(b == 0) return a;
  4     else return gcd(b , a%b);
  5 }
  6 
  7 int lcm(int a , int b)
  8 {
  9     return a/gcd(a , b)*b;//先除后乘防止溢出
 10 }
 11 
 12 inline int my_abs(int x)
 13 {
 14     return x>=0?x:-x;
 15 }
 16 
 17 inline void my_swap(int &a , int &b)
 18 {
 19     int tmp = a;
 20     a = b , b = tmp;
 21 }
 22 /*
 23 高斯消元解线性方程组的解
 24 返回-2表示有浮点数的解，但无整数解
 25 返回-1表示无解
 26 返回0表示唯一解
 27 大于0，表示无穷解，返回的是自由变元个数
 28 有equ个方程，var个变元。增广矩阵行数为equ，分别为0到equ-1，列数为var+1，分别为0到var
 29 总是考虑上三角部分，所以内部更新总是只更新了上三角部分
 30 */
 31 int gauss(int equ , int var)
 32 {
 33     int max_r; // 当前这列绝对值最大的行
 34     int col; // 当前处理的列
 35     int ta , tb , i , j , k;
 36     int Lcm , tmp , free_x_num , free_index;
 37 
 38     for(i=0 ; i<=var ; i++)
 39         x[i] = 0 , free_x[i] = true;
 40 
 41     //转化为阶梯形矩阵
 42     col = 0;//一开始处理到第0列
 43     for(k=0 ; k<equ && col < var; k++,col++){
 44     /*
 45     枚举当前处理的行
 46     找到该行col列元素绝对值最大的那行与第k行交换，
 47     这样进行除法运算的时候就可以避免小数除以大数了
 48     */
 49         max_r = k;
 50         for(i = k+1 ; i<equ ; i++)
 51             if(my_abs(a[i][col]) > my_abs(a[max_r][col]))
 52                 max_r = i;
 53         if(max_r != k){
 54             //交换两行的数据
 55             for(i=k ; i<var+1 ; i++)
 56                 my_swap(a[max_r][i] , a[k][i]);
 57         }
 58 
 59         if(a[k][col] == 0){
 60             //说明该col列第k行以下全是0了，则处理当前行的下一列
 61             k--;
 62             continue;
 63         }
 64 
 65         for(i=k+1 ; i<equ ; i++){
 66             //枚举要删去的行
 67             if(a[i][col] != 0 ){
 68                 Lcm = lcm(my_abs(a[i][col]) , my_abs(a[k][col]));
 69                 ta = Lcm / my_abs(a[i][col]);
 70                 tb = Lcm / my_abs(a[k][col]);
 71                 if(a[i][col] * a[k][col] < 0) tb = -tb;//异号相加
 72                 for(j = col ; j<var+1 ; j++)
 73                     a[i][j] = a[i][j]*ta - a[k][j]*tb ;
 74             }
 75         }
 76     }
 77 
 78     //1.无解的情况
 79     for(i=k ; i<equ ; i++){
 80         if(a[i][col] != 0) return -1;
 81     }
 82 
 83     //2.无穷解：在var*(var+1) 的增广阵中出现了(0,0,0....,0)这样的行，也即说明没有出现严格的上三角
 84     //出现这样情况的行数也就是自由变元的个数
 85     if(k < var){
 86         //首先初始化自由变元有var-k个，即不确定的变元至少有var-k个
 87         for(i=k-1 ; i>=0 ; i--){
 88             //第i行一定不会是(0,0,0,...,0)的情况，因为这样的行实在第k行到equ行
 89             //同样，第 i 行一定不会是(0,0 ... , a) , a != 0的情况
 90 
 91             free_x_num = 0;
 92             //free_x_num记录自由变元个数，如果超过1个，则解集无穷，无法求解.
 93 
 94             for(j=0 ; j<var ; j++){
 95                 if(a[i][j] != 0 && free_x[j]) free_x_num++ , free_index = j;
 96             }
 97 
 98             if(free_x_num > 1) continue; //无法求解出确定变元
 99 
100             tmp = a[i][var];
101             for(j=0 ; j<var ; j++)
102                 if(a[i][j] != 0 && j != free_index)
103                     tmp -= a[i][j]*x[j];
104             x[free_index] = tmp / a[i][free_index];//求出该变元
105             free_x[free_index] = false;//该变元确定
106         }
107         return var - k;//自由变元有var-k个
108     }
109 
110     //唯一解的情况：在var*(var+1) 的增广矩阵中形成严格的上三角阵
111     //计算x[n-1] , x[n-2] .... x[0] 回代计算，先得到后面的值，推到前面的值
112     for(int i=var-1 ; i>=0 ; i--){
113         tmp = a[i][var];
114         for(int j=i+1 ; j<var ; j++){
115             if(a[i][j] != 0 ) tmp -= a[i][j]*x[j];
116         }
117        // if(tmp % a[i][i] != 0) return -2; //说明浮点数存在，无整数解，但是有小数解
118         x[i] = tmp / a[i][i];
119     }
120     return 0;
121 }

 

 posted on 2015-01-21 15:26  Love风吟  阅读(264)  评论(0)  编辑  收藏  举报