Solve a set of simultaneous linear equations by Gauss Elimination with maximization of pivot element(線上執行) (原始碼下載)
對於解決一般的線性聯立方程組，傳統的 Gauss Elimination 或 Gauss-Jordan Elimination 就可應付，這部分上次已經介紹過了；但是如果遇到主對角線元素比其它元素小很多，或是主對角線元素等於零的線性聯立方程組，就必須用到 maximization of pivot elements 的技巧，因為主對角線元素很小將會產生嚴重的捨入誤差，導致計算結果錯誤的問題；主對角線元素為零，將會產生溢位的問題。
maximization of pivot elements 的演算法大致有三種 maximal column pivoting、scaled-column pivoting 以及 complete pivoting，以下要介紹的是 complete pivoting 應用於 Gauss Elimination 的演算法，complete pivoting 是三種方法中執行效率最差，但卻是可將捨入誤差降至最小的方法。(其實我對於 pivoting 並不很熟悉，如果說明有誤，請來信指正，謝謝！)
complete pivoting 演算法是將係數矩陣中絕對值最大的元素找出，稱此元素為 pivot element，一個 n×n 的係數矩陣必須找出 n 個 pivot element，每個 pivot element 不可同列或同行，因為若有兩個pivot element 同列或同行將無法進行 Gauss Elimination，以下以一個簡單的線性聯立方程組來說明 complete pivoting 應用於 Gauss Elimination 的演算法。
step1：在係數矩陣中絕對值最大的元素為 9，所以將第三列與第一列對調，以便將 9 移到(1,1)的位置。
step2：將第一行與第二行對調，以便將 9 移到(1,1)的位置。
step3：將第一列每個元素同除 9。
step4：執行消去法使得第一行的其餘元素為零。
step5：在係數矩陣中絕對值最大的元素為 71/9 (此時不用考慮第一列以及第一行的元素)，所以將第三列與第二列對調，以便將 71/9 移到(2,2)的位置。
step6：將第二行與第三行對調，以便將 71/9 移到(2,2)的位置。
step7：將第二列每個元素同除 71/9。
step8：執行消去法使得第三列第二行的元素為零。
重複以上步驟，直到係數矩陣化簡為上三角矩陣。這裡只是說明如何選取 pivot element，所以並沒有算出最後的結果。
副程式
'以下是Gauss elimination with maximization of pivot的副程式 Public Sub GEMP(m() As Single, ReturnValue() As Single) Dim i As Long, j As Long, k As Long, n As Long Dim r As Long, c() As Long, row As Long, col As Long Dim Pivot As Single, temp() As Single row = UBound(m, 1): col = UBound(m, 2) If col <> row + 1 Then MsgBox "矩陣輸入有誤": Exit Sub ReDim c(1 To row), ReturnValue(1 To row), temp(1 To col) For i = 1 To row Pivot = 0 For j = i To row For k = i To row If Abs(m(k, j)) > Pivot Then Pivot = Abs(m(k, j)) r = k: c(i) = j End If Next k Next j If r <> i Then For j = 1 To col temp(j) = m(i, j) m(i, j) = m(r, j) m(r, j) = temp(j) Next j End If If c(i) <> i Then For j = 1 To row temp(j) = m(j, i) m(j, i) = m(j, c(i)) m(j, c(i)) = temp(j) Next j End If Pivot = m(i, i) If Pivot <> 1 Then For k = 1 To col m(i, k) = m(i, k) / Pivot Next End If For j = i + 1 To row Pivot = m(j, i) If Pivot <> 0 Then For k = 1 To col m(j, k) = m(j, k) - m(i, k) * Pivot Next End If Next Next For i = row To 1 Step -1 ReturnValue(i) = m(i, col) For j = i + 1 To row ReturnValue(i) = ReturnValue(i) - m(i, j) * ReturnValue(j) Next Next For i = row To 1 Step -1 If c(i) <> i Then temp(1) = ReturnValue(c(i)) ReturnValue(c(i)) = ReturnValue(i) ReturnValue(i) = temp(1) End If Next i End Sub
參數說明：
- m：線性聯立方程組之增廣矩陣。
- ReturnValue：傳回未知數的計算結果。
範例一：主對角線元素比其它元素小很多的線性聯立方程組
考慮一個如下的線性聯立方程組
0.0001 10000 100000.0001 3 5 53正解為 x1=1、x2=10，如果沒有使用 maximization of pivot elements 的技巧，將會解出 x1=0、x2=10 (變數宣告為 Single 才會這樣，若變數宣告為 Double 可得正解。) 因為電腦會把 100000.0001 除以 0.0001 算成 1E+09，少了1，因此將得到錯誤的結果。使用 GEMP 副程式將可避免這種錯誤，請參考以下程式：
Private Sub Command1_Click() Dim i As Long, Keyin As String, m() As Single, x() As Single Keyin = InputBox("請輸入線性聯立方程式的增廣矩陣(augmented matrix)，" _ & "以分號"";""來分隔列元素，以空白來分隔行元素。" _ & vbNewLine & "例如有一3×4的矩陣:" & vbNewLine & "1 2 3 4" & vbNewLine & _ "5 6 7 8" & vbNewLine & "9 10 11 12" & vbNewLine & _ "則輸入: 1 2 3 4;5 6 7 8;9 10 11 12") If SepStrToMatrix(Keyin, ";", " ", m) Then Call GEMP(m, x) For i = 1 To UBound(x) Debug.Print x(i) Next Else MsgBox "矩陣輸入有誤，請重新輸入。" End If End Sub其中 SepStrToMatrix 這個副程式，請參閱本站的 將字串拆解成陣列。
執行後在出現輸入盒時輸入 0.0001 10000 100000.0001;3 5 53，最後將在VB的即時運算視窗顯示以下結果：
1 10
範例二：主對角線元素等於零的線性聯立方程組
以下是某位網友問我的問題，這是一個 14×15 的矩陣。(不足的元素就是 0)
1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 3.857 0 0 0 1 -0.516 0 0 0 0 0 0 1 3.857 0 0 0 0 0 0 0 0 0 1 -1.937 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 -1要解出這個線性聯立方程組，如果沒有使用 maximization of pivot elements 的技巧，一定會發生溢位的問題。假設以上矩陣儲存在 C:\test\matrix.txt，執行以下的程式將可在VB的即時運算視窗顯示以下結果：
1 -0.2592689 -0.7407311 0.2592689 0.5024591 -0.761728 -1.000508 0.2594006 0.7411077 -0.5024592 -0.2594007 0.7618599 -2.502967 1.502967Private Sub Command2_Click() Dim i As Long, FileNo As Long, Keyin As String Dim temp As String, m() As Single, x() As Single FileNo = FreeFile Open "c:\test\matrix.txt" For Input As #FileNo Do While Not EOF(FileNo) Line Input #FileNo, temp If Trim(temp) <> "" Then Keyin = Keyin & temp & ";" Loop Close #FileNo Keyin = Mid(Keyin, 1, Len(Keyin) - 1) If SepStrToMatrix(Keyin, ";", " ", m) Then Call GEMP(m, x) For i = 1 To UBound(x) Debug.Print x(i) Next Else MsgBox "矩陣輸入有誤，請重新輸入。" End If End Sub其中 SepStrToMatrix 這個副程式，請參閱本站的 將字串拆解成陣列。
要知道解出的結果是否正確，最笨的方法就是將每個未知數代回方程式中驗算，但是係數矩陣若是 Singular，將可能有無限多組解，因此為了謹慎起見，通常在解線性聯立方程組時會一併計算係數矩陣的 Determinant，以確保該係數矩陣不是 Singular，下次有空再介紹 Determinant 的計算。
GEMP 和 GJEMP 兩個副程式都是用於解決線性聯立方程組的問題，主要差異在於前者是使用 Gauss Elimination，後者是使用 Gauss-Jordan Elimination，由於 GEMP 運算次數少，所以執行效率快。
其實 GEMP 副程式寫的並不好，因為在尋找絕對值最大的元素時是使用最簡單的循序搜尋法，應該有更好的搜尋法可以利用。
This page was written by Jaric on Oct. 10, 1998. All rights reserved.
Total pageview since 4/6/1999.