Visual Basic 數值分析 - Solve a set of simultaneous linear equations by Gauss Elimination

Solve a set of simultaneous linear equations by Gauss Elimination. (線上執行) (原始碼下載)

考慮一個如下的線性聯立方程組：

a_₁₁x_₁+a_₁₂x_₂+...+a_{_1n}x_{_n}=b_₁
a_₂₁x_₁+a_₂₂x_₂+...+a_{_2n}x_{_n}=b_₂
.
.
.
a_{_n1}x_₁+a_{_n2}x_₂+...+a_{_nn}x_{_n}=b_{_n}

由矩陣乘法的定義可知，上述的線性聯立方程組可寫成 AX=B，其中 A 稱為係數矩陣。

       a_₁₁   a_₁₂    ...   a_{_1n}
       a_₂₁   a_₂₂   ...   a_{_2n}
A=   .       .    .        .
        .       .       .     .
        .       .          . .
      a_{_n1}   a_{_n2}     ...   a_{_nn}

       b_₁
       b_₂
B=   .
        .
        .
       b_{_n}

       x_₁
       x_₂
X=   .
        .
        .
       x_{_n}

若以 Guass Elimination 來求解此系統，則通常以一增廣矩陣(augmented matrix)來表示此系統。

a_₁₁   a_₁₂   ...    a_{_1n}   b_₁
a_₂₁   a_₂₂   ...    a_{_2n}   b_₂
.       .    .        .      .
  .       .       .    .      .
.       .          . .      .
a_{_n1}   a_{_n2}    ...    a_{_nn} b_{_n}

並經由一連串的初等列運算，將增廣矩陣中的係數矩陣部分化簡為上三角矩陣，最後以後向代換法(backward substitution)，求出每個未知數。這裡所討論的系統，是方程式的數目和未知數的數目相等，也就是說係數矩陣 A 為一方陣，因此增廣矩陣一定是 n×(n+1) 之矩陣。在將係數矩陣化簡為上三角矩陣時，若遇到主對角線元素為零，就必須使用初等列運算，使得主對角線元素不為零，否則程式會發生溢位的問題。

副程式：

Public Sub GaussElimination(m() As Single, ReturnValue() As Single)
     Dim i As Long, j As Long, k As Long
     Dim row As Long, col As Long, Pivot As Single
     row = UBound(m, 1)
     col = row + 1
     ReDim ReturnValue(1 To row)
     For i = 1 To row
          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
End Sub

上述的 GaussElimination 副程式至少還有三個問題存在，但是對於解決一般的線性聯立方程組，仍然可以適用。

在副程式中並沒有判斷增廣矩陣 m 是否為 n×(n+1)。
沒有使用 Pivoting strategy，因此無法處理主對角線元素可能為零，以及捨入誤差(roundoff error)很嚴重的問題。關於包含 Pivoting strategy 的 GaussElimination 副程式將在以後介紹。
無法處理 ill-conditioned 的問題。

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

參數說明：

m：線性聯立方程組之增廣矩陣。
ReturnValue：傳回未知數的計算結果。

範例：

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

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

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

今有一線性聯立方程組，其增廣矩陣如下：

       3      -2       7      15
      -2       4      -3      12
      -1       9       4      27

執行 Command1 之後，在 VB5 的即時運算視窗可得此線性聯立方程組的三個未知數：

100.3333
28.66667
-32.66667

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

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

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

若有一線性聯立方程組，其增廣矩陣儲存在"C:\m.txt"

       2      -1       4      -1       5       1       2      24
       1      -1       4       3       7       1       5      53
      -1       2      -1      -7       1      -1       1     -10
       1      12       1      -8       4      -3       7      91
       8       4       5      -4       1      -1       5      27
       2       3       9      -1       4       2       1      43
       4       3       1       2      -1      12      -1      47

執行 Command2 之後，在 VB5 的即時運算視窗可得此線性聯立方程組的七個未知數。

2.694457
9.022199
-2.15419
4.831712
8.757999
0.7240171
-1.71613

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

Gauss Elimination 和 Gauss-Jordan Elimination 同樣可以解決線性聯立方程組的問題，但是 Gauss Elimination 執行速度較快，因為 Gauss-Jordan Elimination 的總運算次數較多。

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

Revised : Oct. 10 , 1998 .

Total pageview since 4/6/1999.