文章目錄
- 概念
- 用法
- 特殊情況
- 我的奇怪方法
概念
什么是高斯消元?讓我們看一看 OI-Wiki 的解釋:
高斯消元法(Gauss–Jordan elimination)是求解線性方程組的經典算法,它在當代數學中有著重要的地位和價值,是線性代數課程教學的重要組成部分。
高斯消元法除了用于線性方程組求解外,還可以用于行列式計算、求矩陣的逆,以及其他計算機和工程方面。
通俗來講,就是解多元一次方程。
用法
高斯消元法的核心思想和它的名字一樣:消元。為了方便講解,現在讓我們寫出一個三元一次方程組:
{ x + y + z = 6 1 ? 5 x + 3 y + 2 z = 17 2 ? 3 x + 2 y + z = 10 3 ? \begin{cases}x+y+z=6\ \textcircled{1}\\5x+3y+2z=17\ \textcircled{2}\\3x+2y+z=10\ \textcircled{3}\end{cases} ? ? ??x+y+z=6?1?5x+3y+2z=17?2?3x+2y+z=10?3??
相信各位讀者一定一眼就看得出答案,但計算機不是人,不像各位讀者那么聰明(就我一個蒟蒻,嗚嗚嗚),于是請各位讀者跟隨我穿越時光,回到初二的課堂。
初中老師講過二元一次方程組的求解方法有兩種:代入消元法和加減消元法。代入消元就是指用一個未知數表示出另一個未知數,然后再代入到另一個方程里,而加減消元法則是通過兩個方程括倍再通過加減法消掉一個未知數。
相信各位讀者已經想到高斯消元到底是什么了。好的,本期講解完畢,我們下次再見!
開個玩笑,我們繼續。
剛才我們說到二元一次方程組的求解方式有加減消元和代入消元,那我們現在把它們擴展到 n n n 元一次方程組,結果會怎么樣呢?
答案:代入消元太復雜,加減消元得高消。
這里的高消就是指高斯消元。
很明顯,代入消元對計算機和程序員來講難度太大了,因此,一位偉大的數學家——高斯,便發明了高斯消元(不知他的初衷是啥)!
高斯消元的做法是這樣的:首先把每一項的系數和右邊的常數項提取出來構成矩陣,比如上面的那個方程就可以被轉化成這樣:
{ x + y + z = 6 5 x + 3 y + 2 z = 17 3 x + 2 y + z = 10 → [ 1 1 1 6 5 3 2 17 3 2 1 10 ] \begin{cases}x+y+z=6\\5x+3y+2z=17\\3x+2y+z=10\end{cases}\to\begin{bmatrix}1&1&1&6\\5&3&2&17\\3&2&1&10\end{bmatrix} ? ? ??x+y+z=65x+3y+2z=173x+2y+z=10?→ ?153?132?121?61710? ?
接著讓我們把第 i i i 行第 i i i 列的數變成 1 1 1,也就是給第 i i i 行的每個數除以 a i , i a_{i,i} ai,i?,在上面的矩陣中我們先處理第一行:
[ 1 1 1 6 5 3 2 17 3 2 1 10 ] → [ 1 1 1 6 5 3 2 17 3 2 1 10 ] \begin{bmatrix}1&1&1&6\\5&3&2&17\\3&2&1&10\end{bmatrix}\to\begin{bmatrix}1&1&1&6\\5&3&2&17\\3&2&1&10\end{bmatrix} ?153?132?121?61710? ?→ ?153?132?121?61710? ?
然后我們以第 i i i 行第 i i i 列的數為基準,把第 i i i 列第 i i i 行下面的所有數變成 0 0 0,這一步解釋起來有點復雜,大家自己看圖:
[ 1 1 1 6 5 3 2 17 3 2 1 10 ] ? [ 1 1 1 6 5 ? 5 × 1 3 ? 5 × 1 2 ? 5 × 1 17 ? 5 × 6 3 ? 3 × 1 2 ? 3 × 1 1 ? 3 × 1 10 ? 3 × 6 ] ? [ 1 1 1 6 0 ? 2 ? 3 ? 13 0 ? 1 ? 2 ? 8 ] \begin{aligned}&\begin{bmatrix}1&1&1&6\\5&3&2&17\\3&2&1&10\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\5-5\times1&3-5\times1&2-5\times1&17-5\times6\\3-3\times1&2-3\times1&1-3\times1&10-3\times6\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&-2&-3&-13\\0&-1&-2&-8\end{bmatrix}\end{aligned} ? ?153?132?121?61710? ?? ?15?5×13?3×1?13?5×12?3×1?12?5×11?3×1?617?5×610?3×6? ?? ?100?1?2?1?1?3?2?6?13?8? ??
最后重復上述操作,直到形成一個上三角結構,為了方便讀者理解,這里把全過程展示了出來:
[ 1 1 1 6 0 ? 2 ? 3 ? 13 0 ? 1 ? 2 ? 8 ] ? [ 1 1 1 6 0 ? 2 ÷ ( ? 2 ) ? 3 ÷ ( ? 2 ) ? 13 ÷ ( ? 2 ) 0 ? 1 ? 2 ? 8 ] ? [ 1 1 1 6 0 1 3 2 13 2 0 ? 1 ? 2 ? 8 ] ? [ 1 1 1 6 0 1 3 2 13 2 0 ? 1 ? 1 × ( ? 1 ) ? 2 ? 3 2 × ( ? 1 ) ? 8 ? 13 2 × ( ? 1 ) ] ? [ 1 1 1 6 0 1 3 2 13 2 0 0 ? 1 2 ? 3 2 ] ? [ 1 1 1 6 0 1 3 2 13 2 0 0 ? 1 2 ÷ ( ? 1 2 ) ? 3 2 ÷ ( ? 1 2 ) ] ? [ 1 1 1 6 0 1 3 2 13 2 0 0 1 3 ] \begin{aligned}&\begin{bmatrix}1&1&1&6\\0&-2&-3&-13\\0&-1&-2&-8\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&-2\div(-2)&-3\div(-2)&-13\div(-2)\\0&-1&-2&-8\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&-1&-2&-8\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&-1-1\times(-1)&-2-\frac{3}{2}\times(-1)&-8-\frac{13}{2}\times(-1)\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&0&-\frac{1}{2}&-\frac{3}{2}\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&0&-\frac{1}{2}\div(-\frac{1}{2})&-\frac{3}{2}\div(-\frac{1}{2})\end{bmatrix}\\&\Rightarrow\begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&0&1&3\end{bmatrix}\end{aligned} ? ?100?1?2?1?1?3?2?6?13?8? ?? ?100?1?2÷(?2)?1?1?3÷(?2)?2?6?13÷(?2)?8? ?? ?100?11?1?123??2?6213??8? ?? ?100?11?1?1×(?1)?123??2?23?×(?1)?6213??8?213?×(?1)? ?? ?100?110?123??21??6213??23?? ?? ?100?110?123??21?÷(?21?)?6213??23?÷(?21?)? ?? ?100?110?123?1?6213?3? ??
于是我們就得到了一個上三角結構的矩陣(即下半部分的三角形全是 0 0 0),而它的等效方程就是這個:
[ 1 1 1 6 0 1 3 2 13 2 0 0 1 3 ] ? { x + y + z = 6 y + 3 2 z = 13 2 z = 3 \begin{bmatrix}1&1&1&6\\0&1&\frac{3}{2}&\frac{13}{2}\\0&0&1&3\end{bmatrix}\Leftrightarrow\begin{cases}x+y+z=6\\y+\frac{3}{2}z=\frac{13}{2}\\z=3\end{cases} ?100?110?123?1?6213?3? ??? ? ??x+y+z=6y+23?z=213?z=3?
接下來就是最重要的一步了:代入消元。有人會問:坐著你前面不是說代入消元很麻煩嗎,怎么現在又要代入消元了?是不是玩不起啊?
其實不是的,因為通過觀察左邊的那個方程,我們會發現 z z z 的值實際上已經解開了,對應到矩陣上就是最后一行的倒數第二個是一個 1 1 1,而前面的數全是 0 0 0,這時最后一行的最后一個數就是最后的那個未知數的解,記錄下來。
然后回到倒數第二行,我們不難發現:除了倒數第二個數(注意:我這里沒寫第二個,而寫的是倒數第二個)是 1 1 1 之外,其他的不是 0 0 0 就是其他亂七八糟的數,在右邊的方程中,我們在把 z z z 解開后肯定是帶入到倒數第二個方程中,解開 y y y,而對應在左邊的矩陣里就是用最后一個數減去最后那個數乘以下面那個未知數的值,于是得到了 y = 2 y=2 y=2 這個解,記錄下來。
然后就不需要我講了吧。(不會有人不會舉一反三吧。)
特殊情況
上面的操作中,由于我們過分理想化的假設,導致我們沒遇到一種特殊情況:系數為 0 0 0!我們舉個例子看看會發生什么:
{ y + z = 3 x + z = 4 x + y = 5 \begin{cases}y+z=3\\x+z=4\\x+y=5\end{cases} ? ? ??y+z=3x+z=4x+y=5?
畫成矩陣:
[ 0 1 1 3 1 0 1 4 1 1 0 5 ] \begin{bmatrix}0&1&1&3\\1&0&1&4\\1&1&0&5\end{bmatrix} ?011?101?110?345? ?
然后按照上面的步驟:首先我們需要將 a 1 , 1 a_{1,1} a1,1? 變成 1 1 1。但我們發現個事: a 1 , 1 = 0 a_{1,1}=0 a1,1?=0!按照正常的做法只需要除以一個 a 1 , 1 a_{1,1} a1,1? 就行了,但是沒有除以 0 0 0 這個東西啊喂!于是你的程序會直接 R E \color{purple}{RE} RE 掉。
但是這個方程沒有解嗎?并不是。上過小學奧數或者初中的人都知道這個方程有解 { x = 3 y = 2 z = 1 \begin{cases}x=3\\y=2\\z=1\end{cases} ? ? ??x=3y=2z=1?(呀,一不小心說出答案了),那這是怎么一回事呢?
原因很簡單:只是因為這個方程剛好系數為 0 0 0。你看下一個方程系數不就又是 1 1 1 了嗎?因此我們可以適當把第二行和第一行的數換個位置,于是這個矩陣變成了這樣:
[ 1 0 1 4 0 1 1 3 1 1 0 5 ] \begin{bmatrix}1&0&1&4\\0&1&1&3\\1&1&0&5\end{bmatrix} ?101?011?110?435? ?
然后就可以快樂的解方程了。
總結一下: 規則就兩條:
- 能直接解的直接解。
- 不能直接解的換兩行再解。
最后說明一下判無解和無數解的辦法:無解就是存在某一行除了最后一列其他全是 0 0 0,無數解就是存在一行全是 0 0 0。
淺淺算一下時間復雜度: O ( n 3 ) O(n^3) O(n3)。
代碼太簡單,就不展示了。(我絕對不會告訴你是因為我不會寫。)
但是!
我的奇怪方法
由于本蒟蒻太菜,實在不會寫高消,于是我自己再高消的基礎上改了一下。
首先,我不需要換行,對于每一行我都遍歷一遍找一個主元,而不是一定以第 i i i 行的第 i i i 個為主元。
其次,我每次消是把這一列全消成 0 0 0(除這個主元之外),這樣最后就是一個 01 01 01 矩陣。
最后,我不需要還原,因為左邊是個 01 01 01 矩陣,所以系數只有可能是 0 0 0 或 1 1 1,所以直接輸出最后一列就行。
當然,我的解是隨機分布的,因此需要拿個標記數組標記每個解再第幾行再輸出。
時間復雜度: O ( n 3 ) O(n^3) O(n3)。
代碼:
//洛谷P3389
#include<bits/stdc++.h>
#define int long long
#define code using
#define by namespace
#define plh std
code by plh;
int n,vis[106],h[106];
const double eps=1e-7;
double a[106][106],f[106];
void Gauss()
{for(int i=1;i<=n;i++){int id=0;for(int j=1;j<=n;j++){if(fabs(a[i][j])>eps&&!vis[j]){id=j;break;}}if(id){vis[id]=1;h[id]=i;for(int j=1;j<=n;j++){if(j!=id){a[i][j]/=a[i][id];}}f[i]/=a[i][id];a[i][id]=1;for(int j=1;j<=n;j++){if(j!=i){for(int k=1;k<=n;k++){if(k!=id){a[j][k]-=a[i][k]*a[j][id];}}f[j]-=f[i]*a[j][id];a[j][id]=0;}}}else{if(fabs(f[i])>eps){cout<<"No Solution";n=0;return;}}}
}
signed main()
{cin>>n;for(int i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}cin>>f[i];}Gauss();if(!n){return 0;}for(int i=1;i<=n;i++){if(!vis[i]){cout<<"No Solution";return 0;}}for(int i=1;i<=n;i++){printf("%0.2lf\n",f[h[i]]);}return 0;
}
說句實話,這篇文章制作非常不易(特別是 LateX 那一部分,真的是一個字一個字打出來的),希望讀者能點個贊,非常感謝!