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
0 件のコメント:
コメントを投稿