Visual Basic 數值分析 - Matrix Inversion by Gauss-Jordan Elimination

Matrix Inversion by Gauss-Jordan Elimination (線上執行) (原始碼下載)

一個 n×n 的矩陣 A，若存在一個 n×n 的矩陣 A^{^-1}，使得 AA^{^-1}=A^{^-1}A=單位矩陣，則稱 A^{^-1} 為 A 的逆矩陣，A 稱為非奇異的(nonsingular)，逆矩陣具有唯一性。一個不具逆矩陣的矩陣稱為奇異的(singular)。要求出一個矩陣的逆矩陣有許多方法，這裡介紹以 Gauss-Jordan Elimination 來求出逆矩陣的方法，其他還有一些執行效率比較好的方法(例如：Exchange method)，以後有機會再介紹。

副程式：

Public Sub MatrixInversion(m() As Single, ReturnValue() As Single)
     Dim i As Long, j As Long, k As Long, row As Long, Pivot As Single
     row = UBound(m, 1)
     ReDim ReturnValue(1 To row, 1 To row)
     For i = 1 To row
          ReturnValue(i, i) = 1
     Next
     For i = 1 To row
          Pivot = m(i, i)
          If Pivot <> 1 Then
               For k = 1 To row
                    m(i, k) = m(i, k) / Pivot
                    ReturnValue(i, k) = ReturnValue(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 row
                         m(j, k) = m(j, k) - m(i, k) * Pivot
                         ReturnValue(j, k) = ReturnValue(j, k) - ReturnValue(i, k) * Pivot
                    Next
               End If
          Next
     Next
End Sub

上述的 MatrixInversion 副程式至少還有三個問題存在，但是對於解決一般的矩陣求逆問題，仍然可以適用。

在副程式中並沒有判斷矩陣 m 是否為 n×n (方陣)。
沒有使用 Pivoting strategy，因此無法處理主對角線元素可能為零，以及捨入誤差(roundoff error)很嚴重的問題。
如果 m 為一奇異矩陣，程式將會發生溢位的問題。

P.S. 以上的問題在線上執行中都解決了！

參數說明：

m：欲求逆矩陣之矩陣。
ReturnValue：傳回計算後之逆矩陣。

範例：

1.由鍵盤輸入矩陣元素：

Private Sub Command1_Click()
     Dim Keyin As String, m() As Single, x() As Single
     Keyin = InputBox("請輸入方陣，以分號"";""來分隔列元素，以空白來分隔欄元素。" _
               & vbNewLine & "例如有一3×3的方陣:" & vbNewLine & "1 2 3 " & vbNewLine & _
               "4 5 6 " & vbNewLine & "7 8 9 " & vbNewLine & "則輸入: 1 2 3;4 5 6;7 8 9 ")
     If SepStrToMatrix(Keyin, ";", " ", m) Then
          Call MatrixInversion(m, x)
          Call PrintMatrix(x)
     Else
          MsgBox "方陣輸入有誤，請重新輸入。"
     End If
End Sub

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

若有一矩陣如下：

       1       3       0       1
       2       2       1       3
       3       4       5       6
       0       7       3       1

執行 Command1 之後，在 VB5 的即時運算視窗可得此矩陣之逆矩陣。

7.125 -7.75  3.125 -2.625 
-0.25  0.5 -0.25  0.25 
 2.375 -3.25  1.375 -0.8750001 
-5.375  6.25 -2.375  1.875

2.由檔案讀入矩陣元素：

Private Sub Command2_Click()
     Dim FileNo As Long, Keyin As String, temp As String, m() As Single, x() As Single
     FileNo = FreeFile
     Open "c:\mi.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 MatrixInversion(m, x)
          Call PrintMatrix(x)
     Else
          MsgBox "方陣輸入有誤，請重新輸入。"
     End If
End Sub

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

若有一矩陣儲存在"C:\mi.txt"

       1       2      -1      
       2       1       0  
      -1       1       2

執行 Command2 之後，在 VB5 的即時運算視窗可得此矩陣之逆矩陣。

-0.2222222  0.5555556 -0.1111111 
 0.4444444 -0.1111111  0.2222222 
-0.3333333  0.3333333  0.3333333

MatrixInversion 副程式中的 " If Pivot <> 1 Then "，是為了略去除數為 1 的除法運算，還有 " If i<>j And Pivot <> 0 Then "，是為了略去乘數為 0 的乘法運算，最終目的都是為了提昇執行速度。

如果是奇異矩陣求逆之問題，使用 MatrixInversion 副程式將會出現溢位之問題，或是算出數值很大之逆矩陣，但是這個逆矩陣是錯誤的結果，你可以將原矩陣和算出之逆矩陣相乘，看看是否為單位矩陣，以確保計算結果之正確性，這點要非常注意。

另外，AX=B 之線性方程組也可以用逆矩陣的方式解出 X，X=A^{^-1}B，但是並不建議使用這種方式做，因為逆矩陣的計算效率比較差。

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

Revised : Oct. 10 , 1998 .

Total pageview since 4/6/1999.