高斯消元法


自學了一陣高斯消元啦,感覺這個東西聽着高深,其實還是很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還講了個在線求高斯,大家可以去看看他的博客,概率什么的,查查題目就能找到啦(有趣的“驅趕豬玀”....))

 

好啦....在下就為大家陳述到此。還不懂的歡迎評論~我一定傾盡所能幫助大家。

 

最后附一個解方程整數解的丑代碼....

 

/* 高斯消元法(解整數解 & 化成上三角矩陣來求) */ #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[i][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[i][k],matrix[i][k]=matrix[j][k],matrix[j][k]=t; } int gauss(){ for(int i=1;i<=b_l_cnt;i++) rec_x[i]=0,free_x[i]=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[i]); } else if(flag>0){ printf("It has %d free x\n",flag); for(int i=1;i<=b_l_cnt;i++)                //輸出已知的解 & 還不知道的解 
            if(free_x[i]) printf("x%d is unsure!\n",i); else printf("x%d=%d\n",i,rec_x[i]); } else printf("NO ANSWER"); } int main(){ prework(); mainwork(); return 0; }
View Code

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM