Visual Basic 數值分析 - Solve a set of simultaneous linear equations by Gauss-Jordan Elimination with maximization of pivot elements

Solve a set of simultaneous linear equations by Gauss-Jordan Elimination with maximization of pivot elements. (線上執行) (原始碼下載)

對於解決一般的線性聯立方程組，傳統的 Gauss Elimination 或 Gauss-Jordan Elimination 就可應付，這部分上次已經介紹過了；但是如果遇到主對角線元素比其它元素小很多，或是主對角線元素等於零的線性聯立方程組，就必須用到 maximization of pivot elements 的技巧，因為主對角線元素很小將會產生嚴重的捨入誤差，導致計算結果錯誤的問題；主對角線元素為零，將會產生溢位的問題。

maximization of pivot elements 的演算法大致有三種 maximal column pivoting、scaled-column pivoting 以及 complete pivoting，以下要介紹的是 complete pivoting 應用於 Gauss-Jordan Elimination 的演算法，complete pivoting 是三種方法中執行效率最差，但卻是可將捨入誤差降至最小的方法。(其實我對於 pivoting 並不很熟悉，如果說明有誤，請來信指正，謝謝！)

complete pivoting 演算法是將係數矩陣中絕對值最大的元素找出，稱此元素為 pivot element，一個 n×n 的係數矩陣必須找出 n 個 pivot element，每個 pivot element 不可同列或同行，因為若有兩個pivot element 同列或同行將無法進行 Gauss-Jordan Elimination，以下以一個簡單的線性聯立方程組來說明 complete pivoting 應用於 Gauss-Jordan Elimination 的演算法。

	3 -2 7 15 -2 4 -3 12 -1 9 4 27	step1 →	-1 9 4 27 -2 4 -3 12 3 -2 7 15
step2 →	9 -1 4 27 4 -2 -3 12 -2 3 7 15	step3 →	1 -1/9 4/9 3 4 -2 -3 12 -2 3 7 15
step4 →	1 -1/9 4/9 3 0 -14/9 -43/9 0 0 25/9 71/9 21	step5 →	1 -1/9 4/9 3 0 25/9 71/9 21 0 -14/9 -43/9 0
step6 →	1 4/9 -1/9 3 0 71/9 25/9 21 0 -43/9 -14/9 0	step7 →	1 4/9 -1/9 3 0 1 25/71 189/71 0 -43/9 -14/9 0

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。

重複以上步驟，直到係數矩陣化簡為單位矩陣。這裡只是說明如何選取 pivot element，所以並沒有算出最後的結果。

副程式

'以下是Gauss-Jordan elimination with maximization of pivot的副程式
Public Sub GJEMP(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 = 1 To row
               Pivot = m(j, i)
               If i <> j And Pivot <> 0 Then
                    For k = 1 To col
                         m(j, k) = m(j, k) - m(i, k) * Pivot
                    Next
               End If
          Next
     Next i
     For i = 1 To row
          ReturnValue(i) = m(i, col)
     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

正解為 x_₁=1、x_₂=10，如果沒有使用 maximization of pivot elements 的技巧，將會解出 x_₁=0、x_₂=10 (變數宣告為 Single 才會這樣，若變數宣告為 Double 可得正解。) 因為電腦會把 100000.0001 除以 0.0001 算成 1E+09，少了1，因此將得到錯誤的結果。使用 GJEMP 副程式將可避免這種錯誤，請參考以下程式：

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 GJEMP(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的即時運算視窗顯示以下結果：

Private 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 GJEMP(m, x)
          For i = 1 To UBound(x)
               Debug.Print x(i)
          Next
     Else
          MsgBox "矩陣輸入有誤，請重新輸入。"
     End If
End Sub

其中 SepStrToMatrix 這個副程式，請參閱本站的將字串拆解成陣列。

要知道解出的結果是否正確，最笨的方法就是將每個未知數代回方程式中驗算，但是係數矩陣若是 Singular，將可能有無限多組解，因此為了謹慎起見，通常在解線性聯立方程組時會一併計算係數矩陣的 Determinant，以確保該係數矩陣不是 Singular，下次有空再介紹 Determinant 的計算。

GJEMP 和 GEMP 兩個副程式都是用於解決線性聯立方程組的問題，主要差異在於前者是使用 Gauss-Jordan Elimination，後者是使用 Gauss Elimination，由於 GEMP 運算次數少，所以執行效率快。

其實 GJEMP 副程式寫的並不好，因為在尋找絕對值最大的元素時是使用最簡單的循序搜尋法，應該有更好的搜尋法可以利用。

[ 上一個 | 首頁 | 數值分析 |下一個 ]

Revised : Oct. 10 , 1998 .

Total pageview since 4/6/1999.