數論小白都能看懂的線性方程組及其解法(高斯消元)

此文章依 CC 4.0 BY-SA 版權協議轉載自 ShineEternal 的博客

-1. 序言

說到線性方程組，大家第一反應大概就是高斯消元，本文將對其詳細講解並配合例題與相關方法為您呈現。

本文因圖文並茂有較多配圖且講解詳細較多，再過多的放置代碼會引起文章的冗長以及閱讀的不適，故只將模板放在上面。

0、由來

遇到形如：

\[\left\{\begin{array}{l}{ a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n}=b_{1}} \\ {a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n}=b_{2}} \\ {\cdots} \\ {\cdots} \\ {\cdots} \\ {a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n}=b_{n}} \end{array}\right. \]

這樣的方程組，你會怎么解呢？

如果是在數學領域中，可能這樣的方程組中方程的個數是屈指可數的，這樣就特別簡單，一通亂搞就出來了。

但是如果在 OI 領域遇到這類問題呢？我們發現，計算機中必須由固定的算法，才能實現某一問題，而人腦能進行綜合性的思考，但是卻無法將大腦中的神經元活動表述到程序中，這是人腦與電腦的區別，也就由此誕生了各種解線性方程組的方法。

說到這里，我想起了一本著名的書：計算機與人腦。馮諾依曼寫的，有興趣者可以自行購買閱讀也能普及一些歷史

這種線性的方程組在工程中也運用廣泛，這種模型中：

• 未知量較多
• 方程個數不同

而在考慮解方程時，要考慮：

• 是否有解？
• 如果有解，是否唯一？
• 如果不唯一，解有什么規律與結構？

由此，線性方程組的解決方法就顯得尤為重要。

下面就來介紹一下兩種常用的方法及相關例題

1、高斯消元

插播一些歷史知識：

該方法以數學家高斯命名，由拉布扎比。伊丁特改進，發表於法國但最早出現於中國古籍《九章算術》，成書於約公元前 \(150\) 年。

高斯：

德國著名數學家、物理學家、天文學家、幾何學家，大地測量學家。享有「數學王子」的美譽，真的是多才多藝，著名的「自然數 \(1\sim100\) 的求和問題」就是他 \(9\) 歲時計算出來的，雖然可能現在對我們來說輕而易舉，但是的確體現了他那時的創新精神，為日后的研究打下了堅實的基礎。

這種是大家比較耳熟能詳的解法了。

那么，這種解法的核心思想是什么呢？

我們先來構想一下，什么樣子的線性方程式我們用電腦也能寫出解法？

顯然，我們如果從某一個式子開始，每一個式子都能得到一個未知數的解，而得到了這個未知數的解，又能代入另一個方程，再得到一個未知數的解。

以此類推。

總結一下，就是一個「三角形分割」（自己瞎起的名字）

還是拿上面那個普遍形式：

\[\left\{\begin{array}{l}{ a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n}=b_{1}} \\ {a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n}=b_{2}} \\ {\cdots} \\ {\cdots} \\ {\cdots} \\ {a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n}=b_{n}} \end{array}\right. \]

但是為了觀察方便，我們假設一個三行四列的式子：（真香）

\[\left\{\begin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \\ {x_{1}+4 x_{2}+7 x_{3}=3} \\ {9 x_{1}+3 x_{2}+2 x_{3}=2}\end{array}\right. \]

此時把系數寫到矩陣里：（沒學過矩陣也沒關系，就當成一個表示就行啦）

\[\left[\begin{array}{llll}{1} & {3} & {4} & {5} \\ {1} & {4} & {7} & {3} \\ {9} & {3} & {2} & {2}\end{array}\right] \]

其實就是 高斯消元模板題的樣例

第一步：對左側的系數進行消元

對

\[\left[\begin{array}{lll} {\color{red}1} & {3} & {4} \\ {\color{green}1} & {\color{blue}4} & {7} \\ {\color{green}9} & {\color{orange}3} & {\color{purple}2}\end{array}\right] \]

進行向「三角形分割」的轉換

即目標狀態為

\[\left[\begin{array}{lll}{?} & {?} & {?} \\ {0} & {?} & {?} \\ {0} & {0} & {?}\end{array}\right] \]

此時，我們設標紅的 \(\color{red}1\) 為主元，對照目標狀態可知，標綠的 \(\color{green}1\) 和 \(\color{green}9\) 需要化 \(0\)。

所以我們只需要將第二行整體減去第一行，就實現了將綠 \(\color{green}1\) 化 \(0\)。

那么 \(9\) 如何轉化呢？

此時就對應的整體減去 \(9\) 倍的第一行就行（可能會有更簡便的系數化 \(0\) 法，但是我們要給計算機設計一個統一的算法才能實現）

於是就得到了：

\[\left[\begin{array}{ccc} {1} & {3} & {4} \\ {0} & {\color{blue}1} & {3} \\ {0} & {\color{orange}{-24}} & {-34}\end{array}\right] \]

然后再以藍 \(\color{blue}1\) 為主元，現在我們需要將 \(\color{orange}{-24}\) 化為 0

於是將第三行整體減去 \(−24\) 倍第二行（其實相當於加 \(24\) 倍）

Q：這里不會把之前第一列化好的 \(0\) 給沖掉嗎？
A: 顯然，在上一步已經保證第一列除主元外系數皆為 \(0\)，所以不會出現上述情況

Q: 萬一上面那個綠 \(\color{green}1\) 為 \(0\) 怎么辦？
A: 這就是比較特殊的主元為 \(0\) 的情況，此時因為每一個方程地位都相等，所以只需將主元所在的行與下方主元位置非 \(0\) 的行互換即可。

追加 Q：萬一下面也沒有呢？
A: 那么方程就無解了，這種情況在后面也會提到。

於是就和

\[\left[\begin{array}{lll}{?} & {?} & {?} \\ {0} & {?} & {?} \\ {0} & {0} & {?}\end{array}\right] \]

這個目標狀態對應起來了。

然后就想摩拳擦掌的解了。

但是發現右邊的 \(b\) 還沒有代入矩陣。

於是

第二步：對增廣矩陣消元

所謂增廣矩陣？

增廣矩陣（又稱擴增矩陣）就是在系數矩陣的右邊添上一列，這一列是線性方程組的等號右邊的值。 -----百度百科

其實就是多出 \(b\) 那一列啦

\[\left[\begin{array}{llll}{1} & {3} & {4} & {5} \\ {1} & {4} & {7} & {3} \\ {9} & {3} & {2} & {2}\end{array}\right] \]

如上

增廣矩陣的操作與左側的系數矩陣完全一樣，剛開始單列系數矩陣只是一個由簡入難的過程，加深對高斯消元的理解。限於篇幅原因，就不再贅述。

第三步：回帶

最后經過一系列的消元，就得到了如下矩陣：

\[\left[\begin{array}{cccc}{1} & {3} & {4} & {5} \\ {0} & {1} & {3} & {-2} \\ {0} & {0} & {38} & {-91}\end{array}\right] \]

然后把矩陣再寫回方程組：

\[\left\{\begin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \\ {0 x_{1}+1 x_{2}+3 x_{3}=-2} \\ {0 x_{1}+0 x_{2}+38 x_{3}=-91}\end{array}\right. \]

這樣我們顯而易見地就可以去掉系數為 \(0\) 的字母：

\[\left\{\begin{array}{l}{x_{1}+3 x_{2}+4 x_{3}=5} \\ {1 x_{2}+3 x_{3}=-2} \\ {38 x_{3}=-91}\end{array}\right. \]

現在結果明了了，從下往上推：

依次可以得到

\[x_3​=\frac{-91}{38}​=−2.39 \]

（保留了兩位小數）

然后已知 \(x_3\)，將其代入中間式子中的 \(x_3\)，同理解得

\[x_2=5.18 \]

最后將 \(x_2​,x_3\)​ 都代回最上面的式子，得：

\[x_1=−0.97 \]

於是，對高斯消元的手動模擬就結束了

顯然，這個在數學考試中會答不完題的，但是對於計算機來說，卻是一個通用的方法。

時間復雜度：\(O(n^3)\)

雖然時間復雜度較高，不管怎樣，我們已經找到了解決線性方程組的計算機通用方法，下面就用代碼實現：

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
const double eps=1e-7 ;
double a[105][105],ans[105];
int main()
{
    int n;
	scanf("%d",&n);
    for(int i=1;i<=n;i++)
    { 
        for(int j=1;j<=n+1;j++)//這里千萬別忘了右側那一列了 
        {
            scanf("%lf",&a[i][j]);
        }
    }
    for(int i=1;i<=n;i++)
	{
        int m=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(a[m][i])<fabs(a[j][i]))//fabs 能取浮點數絕對值
                m=j;//找到當前這一列最大的數字，作為主元消元，這樣能最大限度的避免精度誤差 
        if(fabs(a[m][i])<eps)//這里要判精度，其實就是如果該位置的系數為 $0$ 則無解（double 沒法准確處理這種精度 
		{
            printf("No Solution\n");
            return 0;
        }
        if(i!=m)swap(a[i],a[m]);//如果巧了，當前行不是最大行，那得罪了，你換下去吧，讓未來的主元上來，好直接枚舉之后的方程就行了 
        double d=a[i][i];
        for(int j=i;j<=n+1;j++)
            a[i][j]/=d;//將該位置的系數變為 1 
        for(int j=i+1;j<=n;j++)
		{
            d=a[j][i];
            for(int k=i;k<=n+1;k++)
            {
                a[j][k]-=a[i][k]*d;//將其他的方程用兩式相減的方法減去應當減掉的系數的值 
            }
        }
    }
    ans[n]=a[n][n+1];
    for(int i=n-1;i>=1;i--)//回帶，最后一行直接寫出答案，然后其他行還有等待前面處理出來，從后往前推即可 
	{
        ans[i]=a[i][n+1];//這里不難理解 
        for(int j=i+1;j<=n;j++)
        {
            ans[i]-=a[i][j]*ans[j]; 
    	}
    }
    for(int i=1;i<=n;i++)
    { 
        printf("%.2lf\n",ans[i]);
    }
    return 0;
}

這個模板可以在 P3389 【模板】高斯消元法 提交測試（話說貌似我之前已經放過一遍鏈接了

其他的地方應該都不難理解，就是為什么要讓主元是最大的來避免誤差？

這個是從期望上來證明的，有興趣者闊以看看這一篇 日報，應該就能想出來了吧。

大概就是：

設當前要消元的式子系數為 \(q_{i1},q_{i2},q_{i3},\cdots,q_{in}\)

並假定我們主要針對的系數是 \(q_{i1}\)

則有

\[q_{j n}=q_{j n}-\frac{q_{j 1}}{q_{i 1}} \times q_{i w} \]

此時我們求的最大主元就是 \(q_{i1}\)
，它越大，就說明 \(\dfrac{q_{j1}}{q_{i1}}\) ​​的期望越小，這時較小的數給整個結果帶來的影響也小，那么就不容易造成精度誤差了 QAQ

到此，高斯消元就結束了，但是

還有高斯消元的進階版，能避免回帶來求出答案。

2、高斯——約旦消元

說到避免回帶來求出答案，大家應該能想到這種算法的目的，還是拿樣例來說，就是把方程組變成：

\[\left[\begin{array}{llll} {?} & {0} & {0} & {?} \\ {0} & {?} & {0} & {?} \\ {0} & {0} & {?} & {?} \end{array}\right] \]

這樣，每個方程就能獨立的解了。

那么，我們如何得到這種形式呢？

首先，我們依照朴素的高斯消元不難得到：

\[\left[\begin{array}{llll} {?} & {?} & {?} & {?} \\ {0} & {?} & {?} & {?} \\ {0} & {0} & {?} & {?} \end{array}\right] \]

觀察上下兩個矩陣，不難得到，我們在消元時不僅僅消去主元所在式子下方的式子，而對於上方的式子也應當予以處理。

於是就開始了，下方是原圖：

\[\left[\begin{array}{llll}{1} & {3} & {4} & {5} \\ {1} & {4} & {7} & {3} \\ {9} & {3} & {2} & {2}\end{array}\right] \]

第一次與普通消元類似，即得：

\[\left[\begin{array}{rrrr} {1} & {3} & {4} & {5} \\ {0} & {1} & {3} & {-2} \\ {0} & {-24} & {-34} & {-43} \end{array}\right] \]

但是第二次也要將第一行進行消元，得到：

\[\left[\begin{array}{rrrr} {1} & {0} & {-5} & {11} \\ {0} & {1} & {3} & {-2} \\ {0} & {0} & {38} & {-91} \end{array}\right] \]

然后就是最后一步（這一步是比朴素高斯消元要多的）：

\[\left[\begin{array}{rrrr} {1} & {0} & {0} & {-0.97} \\ {0} & {1} & {0} & {5.18} \\ {0} & {0} & {38} & {-91} \end{array}\right] \]

把矩陣寫回方程組

\[\left\{\begin{array}{l}{x_{1}+0 x_{2}+0 x_{3}=-0.97} \\ {0 x_{1}+x_{2}+0 x_{3}=5.18} \\ {0 x_{1}+0 x_{2}+38 x_{3}=-91}\end{array}\right. \]

系數為 \(0\) 去掉:

\[\left\{\begin{array}{l}{x_{1}=-0.97} \\ {x_{2}=5.18} \\ {38 x_{3}=-91}\end{array}\right. \]

這下子結果顯然：

\[\left\{\begin{array}{l}{x_{1}=-0.97} \\ {x_{2}=5.18} \\ {x_{3}=-2.39}\end{array}\right. \]

下面是模板代碼：

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
double a[105][105];
int main()
{
	int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n+1;j++)
        {
            scanf("%lf",&a[i][j]);
        }
    }
    for(int i=1;i<=n;i++)//枚舉列（項） 
    {
        int m=i;
        for(int j=i+1;j<=n;j++)//選出該列最大系數 
        {
            if(fabs(a[j][i])>fabs(a[m][i]))//fabs 是取浮點數的絕對值的函數
            {
                m=j;
            }
        }
        for(int j=1;j<=n+1;j++)//交換
        {
            swap(a[i][j],a[m][j]);
        }
        if(a[i][i]==0)//最大值等於 $0$ 則說明該列都為 $0$，肯定無解 
        {
            printf("No Solution\n");
            return 0;
        }
        for(int j=1;j<=n;j++)//每一項都減去一個數
        {
            if(j!=i)//不是主元那一項 
            {
                double d=a[j][i]/a[i][i];
                for(int k=i+1;k<=n+1;k++)
                {
                    a[j][k]-=a[i][k]*d;
                }
            }
        }
    }
    for(int i=1;i<=n;i++)//最后的結果系數可能不為 $1$，所以記得消去常數 
    {
        printf("%.2lf\n",a[i][n+1]/a[i][i]);
    }
    return 0;
}

顯然，這一份代碼較上一份簡練一些

3、例題

例題 \(1\)

P2455 [SDOI2006] 線性方程組

題意簡敘：

已知 \(n\) 元線性一次方程組，判斷解的情況。

• 無實數解輸出-1
• 無窮多實數解輸出 0
• 有唯一解，則輸出解（小數點后保留兩位小數）。具體格式見樣例

分析：

這道題只是細化了前面模板題對無法得出唯一解情況的判斷，於是順便講一下無解&&多解時的判斷方法：

（我們知道，只要一個未知數無解或多解就可以說明整個方程的情況）

無解

我們最直接的就是判斷最后一行是否是

\[\left[\begin{array}{llll} {?} & {?} & {?} & {?} \\ {0} & {?} & {?} & {?} \\ {0} & {0} & {?} & {?} \end{array}\right] \]

左邊系數為 \(0\) 而右邊系數不為 \(0\)

但在容易疏漏情況，我們考慮到每一行，都會在其前面的消元中將唯一一個未知數留下來，也就都等同於第一行。

於是，每一行都可以按照最后一行的套路模板來判斷：

if(a[i][i]==0&&a[i][n+1]!=0) 無解

多解

與上面對應的，我們很容易想到：

if(a[i][i]==0&&a[i][n+1]==0) 多解

• 要注意先判斷無解，如果是的話直接跳出，然后輪到多解。

因為多解首先建立在有解的基礎上，如果哪一個未知數無解，整個方程自然不存在解，更不必說多解了。

然后其他的就和模板題完全一樣了

例題 2

（這題我沒有找到出處）

題目：【維他命的配方】

題意：

現有若干種維他命，問能否利用這些維他命配制出適合人需求的各種維生素

輸入格式：

第一行：人需補充的維生素種類數 \(V(1\le V\le 25)\)

第二行：\(V\) 個數，第 \(i\) 個數為 \( V_i\)​，表示人體對第 \(i\) 種維生素的需求量 (\(1\le Vi\le 1000\))

第三行：已知的維他命種類 \(G(1\le G\le15)\)

接下來是一個 \(G\times V\) 的整數矩陣，\(A_{ij}\) ​表示第 \(i\) 種維他命中所含的第 \(j\) 中維生素的含量 (\(1\le A_{ij}\le 000\))

輸出格式：

第一行輸出能否配制，能則「YES」, 否則「NO」

第二行：若能配制則輸出 \(G\) 個整數，Gi

Gi​表示第 \(i\) 中維他命所取的數列，多種方案輸出任意一組。不能配制酒無需輸出此行。

輸入樣例：

4
100 200 300 400
4
50 50 50 50
30 100 100 100
20 50 150 250
50 100 150 200

輸出樣例：

YES
1 1 1 0

分析：

首先要弄清維他命與維生素的區別。

然后就應該不難想到高斯消元。

設需要配制第 \(i\) 種維他命數量為 \(x_i(i\le 15)\)

根據題意可列出：

\[\left\{\begin{array}{l}{50 x_{1}+30 x_{2}+20 x_{3}+50 x_{4}=100} \\ {50 x_{1}+100 x_{2}+50 x_{3}+100 x_{4}=200} \\ {50 x_{1}+100 x_{2}+150 x_{3}+150 x_{4}=300} \\ {50 x_{1}+100 x_{2}+250 x_{3}+200 x_{4}=400}\end{array}\right. \]

• 這第一步很好想到，但是行列的順序卻容易弄錯，仔細觀察就會發現需要把 \(V_i\) 補到矩陣的第 \(G+1\) 行，然后每一列為一個方程計算。

列成矩陣的形式：

\[\left[\begin{array}{lllll}{50} & {30} & {20} & {50} & {100} \\ {50} & {100} & {50} & {100} & {200} \\ {50} & {100} & {150} & {150} & {300} \\ {50} & {100} & {250} & {200} & {400}\end{array} \right]\]

一通消元變成：

\[\begin{array}{ccccc}{1} & {2} & {1} & {2} & {4} \\ {0} & {1} & {3 / 7} & {5 / 7} & {10 / 7} \\ {0} & {0} & {1} & {1 / 2} & {1} \\ {0} & {0} & {0} & {0} & {0}\end{array} \]

這時，我們驚奇的發現：最后一行都是 \(0\)

（呀這怎么辦！）

（剛才的知識怎么學的！不就是個多解情況嗎？！）

好了，多不多解不是關鍵，關鍵在於我們要怎么處理這一情況

但這是關鍵，也並不難。

我們只需懶一些將可能多解的未知數設為 \(0\)，來盡可能的減少運算就行。

具體針對這個樣例來說，就是設 \(x_4=0\)，他的數量與答案無關，並不造成影響。

4、高斯消元注意事項：

• 模板不要背錯
• 為了避免精度問題，最好和 eps 比較而不是和 \(0\) 比較（雖然有時這個基本相同）
• 無解多解的判斷情況以及順序
• 遇到題目轉換成高斯消元時不要興奮的轉換出錯或弄反

5、后記

這一篇高斯消元的文章終於完工了，里面的矩陣配圖和線性方程組的 LaTeX 書寫不易，如果這篇文章對您有幫助，那請不吝賜贊。