Excelで一次方程式

　Excelで一次方程式を解くプログラムを作りました。

　最初、LU分解を作ったのですが、使ってみたら、Pivot選択してないから、全然実用にならなくて。で、普通にGaussの消去法を作りました。
　おまけで、反復法のプログラムも作成しています。反復法も使える問題が限られるので、扱いは少し難しいですね。

　なお、このプログラムは画像処理をやるときの伏線です。

　ソースコードは、ご自由にご利用ください。ただし、趣味のプログラムなので、保証はありません。




Option Explicit

Public Const MINIMUM_VALUE As Double = 0.000001
Public Const MAX_ITERATION  As Long = 100

Dim row As Long '結果を表示する行の位置

'LU分解のサンプル
Public Sub matrix_test_LU()
    Dim i As Long
    Dim j As Long
 
    'A x = b
    Dim N As Long
    Dim MatA() As Double
    Dim VecX() As Double
    Dim VecB() As Double
 
    '------------------------------
    'Initialize
    N = 4
 
    ReDim MatA(N * N - 1) As Double
    ReDim VecX(N - 1) As Double
    ReDim VecB(N - 1) As Double
 
    For i = 0 To N - 1
        For j = 0 To N - 1
            MatA(i * N + j) = (1 + i) ^ j
        Next j
     
        VecB(i) = 0
    Next i
 
    VecB(0) = 1
 
    '------------------------------
    'Calculation
    Dim MatLU() As Double
 
    MatrixLU MatA, MatLU, N
    MatrixSolveLU MatA, VecX, VecB, N
 
    '------------------------------
    'Output
    Sheet1.Cells.Clear
 
    row = 1
 
    show_ABX MatA, VecB, VecX, N
 
    Sheet1.Cells(row, 1) = "LU"
    For i = 0 To N - 1
        For j = 0 To i - 1
            Sheet1.Cells(row + i, 2 + j) = MatLU(i * N + j)
        Next j
     
        Sheet1.Cells(row + i, 2 + i) = 1
     
        For j = i To N - 1
            Sheet1.Cells(row + i, 2 + N + j) = MatLU(i * N + j)
        Next j
    Next i
 
    '------------------------------
End Sub

'Gaussの消去法のサンプル
Public Sub matrix_test_Gauss()
    Dim i As Long
    Dim j As Long
 
    'A x = b
    Dim N As Long
    Dim MatA() As Double
    Dim VecX() As Double
    Dim VecB() As Double
 
    '------------------------------
    'Initialize
    N = 5
 
    ReDim MatA(N * N - 1) As Double
    ReDim VecX(N - 1) As Double
    ReDim VecB(N - 1) As Double
 
    For i = 0 To N - 1
        For j = 0 To N - 1
            MatA(i * N + j) = (i + j + 1) Mod N
        Next j
     
        VecB(i) = 0
    Next i
 
    VecB(0) = 1
 
    '------------------------------
    'Calculation
    Dim MatG() As Double
    Dim MatT() As Double
 
    MatrixGauss MatA, MatG, MatT, N
    MatrixSolveGauss MatA, VecX, VecB, N
 
    '------------------------------
    'Output
    Sheet1.Cells.Clear
 
    row = 1
 
    show_ABX MatA, VecB, VecX, N
 
    Sheet1.Cells(row, 1) = "T, A, G"
    For i = 0 To N - 1
        For j = 0 To N - 1
            Sheet1.Cells(row + i, 2 + j) = MatT(i * N + j)
            Sheet1.Cells(row + i, 2 + N + j) = MatA(i * N + j)
            Sheet1.Cells(row + i, 2 + 2 * N + j) = MatG(i * N + j)
        Next j
    Next i
 
    '------------------------------
End Sub

'反復法のサンプル
'参考 <http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdf>
Public Sub matrix_test_iteration()
    Dim i As Long
    Dim j As Long
 
    'A x = b
    Dim N As Long
    Dim MatA() As Double
    Dim VecX() As Double
    Dim VecB() As Double
 
    '------------------------------
    'Initialize
    N = 3
 
    ReDim MatA(N * N - 1) As Double
    ReDim VecX(N - 1) As Double
    ReDim VecB(N - 1) As Double
 
    For i = 0 To N - 1
        For j = 0 To N - 1
            MatA(i * N + j) = 1
        Next j
     
        MatA(i * N + i) = 10
        VecB(i) = 0
    Next i
 
    VecB(0) = 1
 
    '------------------------------
    'Calculation
    MatrixSolveIteration MatA, VecX, VecB, N
 
    '------------------------------
    'Output
    Sheet1.Cells.Clear
 
    row = 1
 
    show_ABX MatA, VecB, VecX, N
 
End Sub

'行列の表示
Public Sub show_ABX(MatA() As Double, VecB() As Double, VecX() As Double, N As Long)
    Dim i As Long
    Dim j As Long
 
    Sheet1.Cells(row, 1) = "A"
    For i = 0 To N - 1
        For j = 0 To N - 1
            Sheet1.Cells(row + i, 2 + j) = MatA(i * N + j)
        Next j
    Next i
 
    row = row + N
 
    Sheet1.Cells(row, 1) = "B"
    For i = 0 To N - 1
        Sheet1.Cells(row + i, 2) = VecB(i)
    Next i
    row = row + N
 
    Sheet1.Cells(row, 1) = "X"
    For i = 0 To N - 1
        Sheet1.Cells(row + i, 2) = VecX(i)
    Next i
 
    row = row + N
End Sub

'LU分解
'A = L * U
Public Sub MatrixLU(MatA() As Double, ByRef MatLU() As Double, N As Long)
    Dim i As Long
    Dim j As Long
    Dim k As Long
 
    ReDim MatLU(N * N - 1) As Double
 
    '------------------------------
    For i = 0 To N - 1
        For j = 0 To N - 1
            MatLU(i * N + j) = MatA(i * N + j)
        Next j
    Next i
 
    '------------------------------
    For i = 0 To N - 1
        For j = 0 To i - 1
            For k = 0 To j - 1
                MatLU(i * N + j) = MatLU(i * N + j) - MatLU(i * N + k) * MatLU(k * N + j)
            Next k
         
            If Math.Abs(MatLU(j * N + j)) < MINIMUM_VALUE Then GoTo Label1
         
            MatLU(i * N + j) = MatLU(i * N + j) / MatLU(j * N + j)
        Next j
     
        For j = i To N - 1
            For k = 0 To i - 1
                MatLU(i * N + j) = MatLU(i * N + j) - MatLU(i * N + k) * MatLU(k * N + j)
            Next k
        Next j
    Next i
 
    '------------------------------
    Exit Sub
 
Label1:
    MsgBox "Divided by Zero"
End Sub

'LU分解
'B = A * X = L * U * X
Public Sub MatrixSolveLU(MatA() As Double, ByRef VecX() As Double, VecB() As Double, N As Long)
    Dim i As Long
    Dim j As Long
 
    Dim MatLU() As Double
    ReDim MatLU(N * N - 1) As Double
    ReDim VecX(N - 1) As Double
 
    MatrixLU MatA, MatLU, N
 
    '------------------------------
    For i = 0 To N - 1
        VecX(i) = VecB(i)
     
        For j = 0 To i - 1
            VecX(i) = VecX(i) - MatLU(i * N + j) * VecX(j)
        Next j
    Next i
 
    '------------------------------
    For i = N - 1 To 0 Step -1
        For j = i + 1 To N - 1
            VecX(i) = VecX(i) - MatLU(i * N + j) * VecX(j)
        Next j
     
        If Math.Abs(MatLU(i * N + i)) < MINIMUM_VALUE Then GoTo Label1
         
        VecX(i) = VecX(i) / MatLU(i * N + i)
    Next i
 
    '------------------------------
    Exit Sub
 
Label1:
    MsgBox "Divided by Zero"
End Sub

'Gaussの消去法
'T * A = G
Public Sub MatrixGauss(MatA() As Double, ByRef MatG() As Double, ByRef MatT() As Double, N As Long)
    Dim i As Long
    Dim j As Long
    Dim k As Long
    Dim temp As Double
 
    ReDim MatG(N * N - 1) As Double
    ReDim MatT(N * N - 1) As Double
 
    '------------------------------
    For i = 0 To N - 1
        For j = 0 To N - 1
            MatG(i * N + j) = MatA(i * N + j)
            MatT(i * N + j) = 0
        Next j
     
        MatT(i * N + i) = 1
    Next i
 
    '------------------------------
    For i = 0 To N - 1
        'find pivot
        k = i
     
        For j = i + 1 To N - 1
            If Math.Abs(MatG(k * N + i)) < Math.Abs(MatG(j * N + i)) Then
                k = j
            End If
        Next j
     
        'row switching
        For j = 0 To N - 1
            temp = MatG(i * N + j)
            MatG(i * N + j) = MatG(k * N + j)
            MatG(k * N + j) = temp
         
            temp = MatT(i * N + j)
            MatT(i * N + j) = MatT(k * N + j)
            MatT(k * N + j) = temp
        Next j
     
        'row multiplication
        If Math.Abs(MatG(i * N + i)) < 0.000001 Then GoTo Label1
     
        temp = 1 / MatG(i * N + i)
     
        For j = 0 To N - 1
            MatG(i * N + j) = temp * MatG(i * N + j)
            MatT(i * N + j) = temp * MatT(i * N + j)
        Next j
     
        'row addition
        For j = i + 1 To N - 1
            temp = MatG(j * N + i)
         
            For k = 0 To N - 1
                MatG(j * N + k) = MatG(j * N + k) - temp * MatG(i * N + k)
                MatT(j * N + k) = MatT(j * N + k) - temp * MatT(i * N + k)
            Next k
        Next j
    Next i
 
    '------------------------------
    Exit Sub
 
Label1:
    MsgBox "Divided by Zero"
End Sub

'Gaussの消去法
'T * B = T * A * X = G * X
Public Sub MatrixSolveGauss(MatA() As Double, ByRef VecX() As Double, VecB() As Double, N As Long)
    Dim i As Long
    Dim j As Long
 
    Dim MatG() As Double
    Dim MatT() As Double
    ReDim MatG(N * N - 1) As Double
    ReDim MatT(N * N - 1) As Double
    ReDim VecX(N - 1) As Double
 
    MatrixGauss MatA, MatG, MatT, N
 
    '------------------------------
    For i = 0 To N - 1
        VecX(i) = 0
     
        For j = 0 To N - 1
            VecX(i) = VecX(i) + MatT(i * N + j) * VecB(j)
        Next j
    Next i
 
    '------------------------------
    For i = N - 1 To 0 Step -1
        For j = i + 1 To N - 1
            VecX(i) = VecX(i) - MatG(i * N + j) * VecX(j)
        Next j
     
        If Math.Abs(MatG(i * N + i)) < 0.000001 Then GoTo Label1
     
        VecX(i) = VecX(i) / MatG(i * N + i)
    Next i
 
    '------------------------------
    Exit Sub
 
Label1:
    MsgBox "Divided by Zero"
End Sub

'反復法
'B = A * X = ( L + D + U) * X
'X = D^-1 * ( B - L * X - U * X)
'Jacobi         : X_i+1 = D^-1 * ( B - L * X_i - U * X_i)
'Gauss-Seidel   : X_i+1 = D^-1 * ( B - L * X_i+1 - U * X_i)
'epsilon = ( B - A * X_i) * ( B - A * X_i) / ( B * B)
'        = ( D * ( X_i+1 - X_i)) * ( D * ( X_i+1 - X_i)) / ( B * B) (for Jacobi)
Public Sub MatrixSolveIteration(MatA() As Double, ByRef VecX() As Double, VecB() As Double, N As Long)
    Dim i As Long
    Dim j As Long
    Dim k As Long
 
    Dim VecY() As Double
 
    ReDim VecX(N - 1) As Double
    ReDim VecY(N - 1) As Double
 
    Dim epsilon As Double
    Dim b2 As Double
 
    '------------------------------
    '収束性の確認
    For i = 0 To N - 1
        epsilon = Math.Abs(MatA(i * N + i))
     
        For j = 0 To i - 1
            epsilon = epsilon - Math.Abs(MatA(i * N + j))
        Next j
     
        For j = i + 1 To N - 1
            epsilon = epsilon - Math.Abs(MatA(i * N + j))
        Next j
     
        If epsilon < 0 Then
            MsgBox "It's not convergent."
            Exit Sub
        End If
    Next i
 
    '------------------------------
    For i = 0 To N - 1
        If Math.Abs(MatA(i * N + i)) < MINIMUM_VALUE Then GoTo Label1
    Next i
 
    '------------------------------
    b2 = 0
 
    For i = 0 To N - 1
        b2 = b2 + VecB(i) * VecB(i)
    Next i
 
    If b2 < MINIMUM_VALUE Then GoTo Label1
 
    '------------------------------
    '反復法の初期値
    For i = 0 To N - 1
        VecX(i) = VecB(i)
    Next i
 
    '------------------------------
    For i = 0 To MAX_ITERATION - 1
        For j = 0 To N - 1
            VecY(j) = VecB(j)
         
            For k = 0 To j - 1
                'Jacobi method
                'VecY(j) = VecY(j) - MatA(j * N + k) * VecX(k)
             
                'Gauss-Seidel method
                VecY(j) = VecY(j) - MatA(j * N + k) * VecY(k)
            Next k
         
            For k = j + 1 To N - 1
                VecY(j) = VecY(j) - MatA(j * N + k) * VecX(k)
            Next k
         
            VecY(j) = VecY(j) / MatA(j * N + j)
        Next j
     
        '------------------------------
        '収束の判定
        epsilon = 0
     
        For j = 0 To N - 1
            epsilon = epsilon + (MatA(j * N + j) * (VecY(j) - VecX(j))) ^ 2
        Next j
     
        epsilon = Math.Sqr(epsilon / b2)
     
        If epsilon < MINIMUM_VALUE Then
            MsgBox "converge : " + CStr(i)
            Exit Sub
        End If
     
        '------------------------------
        For j = 0 To N - 1
            VecX(j) = VecY(j)
        Next j
     
        '------------------------------
    Next i

    '------------------------------
    Exit Sub

Label1:
    MsgBox "Divided by Zero"
End Sub