高斯消元法 附 程序源代码

蜡笔编程

自学了一阵高斯消元啦，感觉这个东西听着高深，其实还是很Logical（有逻辑的）。下面我就分享一下自己对高斯消元的认识啦，希望也可以帮初学者了解这个算法。

首先我们要清楚：高斯消元的目的在于求线性方程组的解。

所以呢，我们先从一个小小的解方程组的例子开始：



伟大的数学天才们....快告诉我233，我把这个方程组AK了。

还在思索的同学们，加油思考~答案就是（2,3,3）

让我们回顾一下初中的老师们怎么教我们解方程的吧~

（因为我发现只能上传10张照片....所以下面就不用写字板截图啦...()里的数字表示下标）

第一步：②*2-①     消去x(1)得 7x(2)+3x(3)=30        .....④

第二步：③*2-①*3    再消一次x(1)得 -x(2)+7x(3)=18    ......⑤

第三步：④+⑤*7    消去x(2)得 52x(3)=156   ......⑥

解之得x(3)=156/52=3

然后代入④或⑤消元得x(2)=3，再将x(2),x(3)的值代入①式，得x(1)=2

是不是觉得很简单呢？....其实这就是高斯消元的思想啦——加减消元 & 代入消元


当然，高斯消元可是要成为'海贼王'的算法（要能解其他一般性的方程），当然要有一定的通用方法来进行这个步骤。

首先我们先运用到一个叫矩阵的东西，我们取出方程组中的各个变量前的系数来成为矩阵的系数，同一个方程的系数放在同一行，同一个变量的系数放在同一列

就像上面的方程可以转成一个3*4的矩阵 A



其中第4列是常数项，其它三列是方程的系数。现在我们再定义一些关于矩阵的基础行变换。

一、交换变换：Ri<->Rj，表示将Ri与Rj的所有元素对应交换

二、倍法变换：Ri=Ri*k，表示将Ri行的所有元素都乘上一个常数k

三、消去变换：Ri=Ri+Rj*k，表示将Ri行的所有元素对应的加上Rj行元素的k倍


wowow~是不是感觉到矩阵的基础行变换和平时加减消元的关系了呢？

矩阵中的元素就是方程中的系数，如此一来，我们可以通过矩阵变换，将方程中的一些变量的系数消为0，然后就可以求出某个变量的值了。

举个例子吧（用上面的矩阵）：R1：  2 3 1 16  ，R2：1 5 2 23

我们想要消掉第一行中的第一个元素，那么就用R1-R2*2（消去变换）得到新的R1 ：0  -7  -3  -30，当然这就相当于我们加减消元啦。


所以我们最后的矩阵渴望得到一个什么样的状态呢？

一般我们会有两个渴望状态：



模式一：【标准型矩阵】

明显的我们可以发现，每行中系数不为0的变量直接用常数项除以该变量的系数就可以得到这个变量的值了。

我们再来讨论一下如何得到这样的矩阵呢？

首先我们从第一行中选出第一列的元素留下，其他所有行都与第一行消去我选择的这列上的系数。如样例：

【2】    3   1  16     //【】表示消去的是哪一列，将其他行上的这一列消为0

【0】  3.5 1.5 15     .....R2=R2-R1*0.5

【0】 -0.5 3.5  9      .....R3=R3-R1*1.5

上面是求实数解的方法，还有求整数解的方法（消元的时候用最小公倍数消去目标系数）

【2】   3   1     16

【0】   7   3    30       因为lcm(1,2)=2，所以R2=R2*(2/1)-R1*(2/2)=R2*2-R1*1

【0】  -1  7    18       因为lcm(3,2)=6，所以R3=R3*(6/3)-R1*(6/2)=R3*2-R1*3

......其他行依次类推，因为当处理第i列的时候，前i-1列都是只有一个系数不为0的元素的，所以加减消元的时候不会影响到（已经决定的元素都是加减一个0*k，没有影响）希望大家自己手推一下这个最小公倍数的下面几步，看看是否得到下面这个矩阵

2   0   0     4

0   7   0    21

0   0  52  156


模式二：【上三角型矩阵】

我们可以将矩阵化成阶梯状，这样如果我们求出了最后一个方程的解，就可以往回代入消元，将已知的变量代入方程，然后一层一层的解出未知数直到第一行得出所有未知数的值。

如何得到一个上三角矩阵呢？

消元方法和标准矩阵一样，都是使用消去变换，让变量系数变成0，然后求解，区别在于每次不需要对每一行都进行消元而只要对当前行下面的元素消元。如样例（用最小公倍数法）：

Step1

【2】  3  1  16

【0】  7  3   30

【0】 -1  7   18

Step2

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。


这便是两种基本的转换方法。都可以解出线性方程组的解。


下面我们再来考虑一些和方程有关的东西。

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还讲了个在线求高斯，大家可以去看看他的博客，概率什么的，查查题目就能找到啦(有趣的“驱赶猪猡”....)）


好啦....在下就为大家陈述到此。还不懂的欢迎评论~我一定倾尽所能帮助大家。


最后附一个解方程整数解的丑代码....

[mw_shl_code=cpp,true]/*
    高斯消元法（解整数解 & 化成上三角矩阵来求）
*/
#include<cstdio>
#include<cstring>
#include<cstdlib>

using namespace std;

#define maxn 1010

int f_c_cnt,b_l_cnt;
int rec_x[maxn];                //记录 x的值的数组
int matrix[maxn][maxn];            //记录系数的矩阵

bool free_x[maxn];                //判断是否为自由元

void prework(){
    freopen("x.in", "r", stdin);
    scanf("%d%d",&f_c_cnt,&b_l_cnt);        //方程数 & 变量数
    for(int i=1;i<=f_c_cnt;i++)
        for(int j=1;j<=b_l_cnt+1;j++)
            scanf("%d",&matrix[j]);
}

int gcd(int a,int b){
    int t;
    while(b) t=b,b=a%b,a=t;
    return a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

void swap(int i,int j){
    int t;
    for(int k=i;k<=b_l_cnt+1;k++)
        t=matrix[k],matrix[k]=matrix[j][k],matrix[j][k]=t;
}

int gauss(){
   
    for(int i=1;i<=b_l_cnt;i++)
        rec_x=0,free_x=true;
   
    int line_0,now_b_l,max_l;                                //line_0 表示这条线以下的全是 0，刚开始会从 1开始向下扫描
   
    for(line_0=now_b_l=1;line_0<=f_c_cnt && now_b_l<=b_l_cnt;line_0++,now_b_l++){
        max_l=line_0;                                         //在当前变量 now_b_l中系数绝对值最大的
        for(int fuc=line_0+1;fuc<=f_c_cnt;fuc++)             //找到变量系数最大的这行
            if(abs(matrix[max_l][now_b_l])<abs(matrix[fuc][now_b_l]))
                max_l=fuc;
        
        if(max_l!=line_0)    swap(line_0,max_l);                  //将绝对值最大的替换当前行，可以防止此行为 0，如果实数运算可以减小误差
        if(!matrix[line_0][now_b_l]) {line_0--;continue;}     //如果这行绝对值最大的也是0，说明这个变量没有系数了，将这个变量过滤
        
        for(int fuc=line_0+1;fuc<=f_c_cnt;fuc++)
            if(matrix[fuc][now_b_l]){
                int LCM=lcm(abs(matrix[fuc][now_b_l]),abs(matrix[line_0][now_b_l]));    //加减消元的时候，用他们的LCM来消元
                int mul=LCM/matrix[fuc][now_b_l];                //这一行应该乘多少
                int div=LCM/matrix[line_0][now_b_l];            //要消元的行要乘多少
                if((long long) matrix[fuc][now_b_l]*matrix[line_0][now_b_l]<0)    div=-div;    //因为之前一直讨论的是绝对值，这里要判断负数
               
                for(int b_l=now_b_l;b_l<=b_l_cnt+1;b_l++)
                    matrix[fuc][b_l]=matrix[fuc][b_l]*mul-matrix[line_0][b_l]*div;    //加减消元
            }
    }
   
            
    for(int fuc=line_0;fuc<=f_c_cnt;fuc++)        //如果所有点都有解，line_0返回值应该是 f_c_cnt+1，如果line_0以下系数就都是 0了
        if(matrix[fuc][b_l_cnt+1])                //如果系数是0而常数项不是0，那么就是无解
            return -1;
        
    if(line_0<f_c_cnt+1){                        //如果不是所有点都能解出来，那么多解的情况便是出现了变元
        for(int fuc=line_0-1;fuc>=1;fuc--){        //在line_0以上的是系数不为 0的方程，考虑才有意义
            int uns_sum=0,uns_num;                //uns_sum 表示 unsure的变量的数量，uns_num表示哪个unsure的是哪个（些）变量
            for(int b_l=1;b_l<=b_l_cnt;b_l++)
                if(matrix[fuc][b_l] && free_x[b_l])
                    uns_sum++,uns_num=b_l;         
            
            if(uns_sum>1)    continue;            //如果一个方程中有超过一个不知道解的变量，那就解不出了
            
            int ans_c=matrix[fuc][b_l_cnt+1];    //否则可以把这个不确定的量解出来，将常数项定为 ans_c
            for(int b_l=1;b_l<=b_l_cnt;b_l++)
                if(matrix[fuc][b_l] && b_l!=uns_num)     //减去其中已知的变量*它们的系数
                    ans_c-=matrix[fuc][b_l]*rec_x[b_l];
            
            rec_x[uns_num]=ans_c/matrix[fuc][uns_num];     //最后的解就是常数项除以要求的变量的系数
            free_x[uns_num]=false;
        }
        
        return b_l_cnt-line_0+1;                //返回变元的个数
    }
        
    for(int fuc=f_c_cnt;fuc>=1;fuc--){            //除了无解和无数解，那就是唯一解了
        int ans_c=matrix[fuc][b_l_cnt+1];        //同上，将解从最后一个变量开始向上代入消元
        for (int b_l=fuc+1;b_l<=b_l_cnt;b_l++)
                ans_c-=matrix[fuc][b_l]*rec_x[b_l];
        rec_x[fuc]=ans_c/matrix[fuc][fuc];
    }
    return 0;
}

void mainwork(){
    int flag=gauss();
    if(flag==0){
        for(int i=1;i<=b_l_cnt;i++)
            printf("x%d=%d\n",i,rec_x);
    }
    else if(flag>0){
        printf("It has %d free x\n",flag);
        for(int i=1;i<=b_l_cnt;i++)                //输出已知的解 & 还不知道的解
            if(free_x)
                printf("x%d is unsure!\n",i);
            else
                printf("x%d=%d\n",i,rec_x);
    }
    else
        printf("NO ANSWER");
}

int main(){
    prework();
    mainwork();
    return 0;
}[/mw_shl_code]