Download presentation
Presentation is loading. Please wait.
1
Gaussian Elimination 東海大學物理系‧數值分析 施奇廷
2
線性聯立方程式之解 如何解這個聯立方程式？ 高斯消去法可直接求精確解
3
(2)式-(1)式×2
4
(3)式-(1)式×3
5
(4)式+(1)式
6
除(1)式外，所有的x1項皆被消掉，依此法再消去(3)(4)式中的x2項以及(4)中的x3項
7
由(4)式就可以很簡單地解出x4，代回(3)解出x3......依此類推即可解出此聯立方程式
8
或是由(4)式開始往上消去x4各項(3)-(4)×13、(2)+(4)×5、(1)-(4) ×3
9
依此類推，消去(1)式中之x2與(2)式中之x3即可
10
解出最後的答案
11
聯立線性方程式之一般表示法
12
將方程式以增廣矩陣（augmented matrix）表示
再利用前述之消去法可將 aij 部分的矩陣對角化，即可解出此聯立方程式之解。
14
Note: 當aii=0 時 由於在消去的過程中，我們會用到 ajk – aik/aii，要是aii = 0 時，該如何處置？
只要找到第 k 列的 aki 不為零，與第 i 列對調，即可得到新的不為零的 aii
15
高斯消去法流程圖 Start Input A(N,N+1) Final Output Output: Singular Matrix
DO I=1,N STOP YES YES A(I,I)=0 I=N? NO NO NO Exist K>I such That A(K,I)≠0? DO J≠I DO K=1,N+1 YES Exchange I-th and K-th Lines A(J,K)=A(J,K)-A(I,K)/A(I,I) 高斯消去法流程圖
16
Scaling of Time 完成上述程式後，試著增加問題的「維度」（未知數與方程式的數目），測試在不同維度（N=10, 50, 100, 300, 500, 1000）下所需的計算時間 注意：augmented matrix 不再用手動輸入，改用公式：
17
有限精度下可能產生的問題 理論上高斯消去法可以精確地解出任何聯立線性方程式之解，然而電腦的精確度有限，在某些情況下可能會發生問題
由於在對第 k 列做第 i 列高斯消去法時必須用到mki=aki/aii，當軸元素（pivot element, 即 aii）為零而造成發散的情形，我們已經以列交換的方式處理了 然而當 aii 雖然不為零，但是非常小時，小到接近電腦精確度時，mki會變得非常大，這時候可能也會有誤差發生
18
精度問題（續） 例如，將前面的例子所有的實數變數宣告為單精度（REAL*4），然後將初始的矩陣改為：
就數學上而言，此方程式之解應與原來一模一樣，然而由於單精度實數有效位數僅有八位，最後一位會被捨去或進位，此時即會造成誤差
19
執行結果
20
解決方式：軸轉策略（Pivoting） 類似碰到aii=0的狀況，找一列aki≠0，將第k列與第i列互換的策略。如果aii很小，找一列aki較大者，將兩列交換即可 要找哪一列k？當然就找所有的aki中最大的那一個，然後與第i列交換
21
範例 考慮二元聯立方程式： 為了方便討論，我們假設電腦精確度只有四位，以下會被四捨五入
因此：m21=a21/a11=1763.6~1764，第一次高斯消去的結果為：
22
範例（續） 我們可以看到，第二個方程式的係數都變得很大，而且因為只能留下四位有效數字，後面的都被捨去了。E2的精確形式應為： y= 精確的y值為1，但由於精度之故所得之近似值為y=1.001，看起來雖然還算好.... 將 y 值代回 E1 求 x 值，得到 x=-10 但是精確解應該是(x, y)=(10,1) 才對！ 會有這麼大的誤差是因為y的0.1%的誤差在E1式中被放大了59.14/0.003=19713倍！
23
解決方法： 由於a11=0.003 << a21，於是將 E1 與 E2 兩列互調，就不會出現軸元素太小的情況，也就不會發生誤差被放大的情形了 新的方程式，同樣在四位精確度下，可以得到正確的解 (x,y)=(10,1)
24
Puzzle: 如果我們再將上面的例子改寫一下，變成：
這組方程式除了將原來的 E1 等號兩邊各乘以10000以外，其他完全一樣。看似沒有軸元素過小的問題（a11=30），不需進行 pivoting，然而很明顯地，使用一般的高斯消去法的結果會得到錯誤的解 (-10,1.001)，為什麼？應該如何解決？
Similar presentations
