Excel (VBA)

Excel VBAに関するフォーラムです。
  • 解決済みのトピックにはコメントできません。
このトピックは解決済みです。
質問

 
(指定なし : 指定なし)
LU分解について
投稿日時: 18/06/06 12:50:57
投稿者: daipon

長文ですがお許しください。
行列AをLU分解し、さらにLの逆行列Mを算定し、L×M=I(単位行列)になるか確認するプログラムを書いてみましたが、0となる部分にわずかに数値が残ります(10^-16のオーダー)。
これは誤差によるものなのか、計算に不備があるのか、教えてもらえないでしょうか。
 
'三角化法による連立方程式の解法LU変換
Sub LU分解()
 '記号の形式を定義
 Dim a(1 To 4, 1 To 4) As Double, b(1 To 4, 1 To 1) As Double
 Dim L() As Double, U() As Double, X() As Double, Y() As Double, Z() As Double, M() As Double, N() As Double
 Dim i As Long, j As Long
 For i = LBound(a, 1) To UBound(a, 1) ' 配列 a の一次元の最少要素数から最大要素数まで
  For j = LBound(a, 2) To UBound(a, 2) ' 配列 a の二次元の最少要素数から最大要素数まで
   a(i, j) = Cells(i + 1, j + 2)
  Next j
 Next i
 For i = LBound(b, 1) To UBound(b, 1) ' 配列 bの一次元の最少要素数から最大要素数まで
  For j = LBound(b, 2) To UBound(b, 2) ' 配列 bの二次元の最少要素数から最大要素数まで
   b(i, j) = Cells(i + 6, j + 2)
  Next j
 Next i
 '動的配列をセット
 ReDim L(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
 ReDim U(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
 'LU分解Call
 Call LU(a, L, U)
 '下三角行列Lの逆行列MCall
 Call LInv(L, M)
 'LM行列の積Call
 X = L
 Y = M
 Call Prod(X, Y, Z)
 'X1行列出力
 For i = 1 To UBound(Z, 1)
  For j = 1 To UBound(Z, 2)
   Cells(i + 41, j + 2) = Z(i, j)
  Next j
 Next i
End Sub
 
'LU分解サブルーチン
Function LU(a() As Double, L() As Double, U() As Double)
 Dim c() As Double
 Dim i As Long, j As Long, k As Long
 'L,Uに分配前の成分を計算(※成分が0,1,aijのものは後で配分する)
 '@c(i,j)を新たに定義し、a(i,j)を初期値として取り込む
 'Ai=1の領域の行列を計算
 'Bi>=j,j>1の領域の行列を計算
 'Ci<j,j>1の領域の行列を計算
 ' 動的配列をセット
 ReDim c(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
 For i = 1 To UBound(a, 1)
  For j = 1 To UBound(a, 2)
   c(i, j) = a(i, j) '@
  If i = 1 Then 'A
   c(i, j) = c(i, j) / c(i, i)
  ElseIf j > 1 And i >= j Then 'B
   For k = 1 To j - 1
    c(i, j) = c(i, j) - c(i, k) * c(k, j)
   Next k
  ElseIf j > 1 And i < j Then 'C
   For k = 1 To i - 1
    c(i, j) = c(i, j) - c(i, k) * c(k, j)
   Next k
   c(i, j) = c(i, j) / c(i, i)
  End If
  Next j
 Next i
 'LU行列に分配
 For i = 1 To UBound(a, 1)
  For j = 1 To UBound(a, 2)
   If i < j Then
    L(i, j) = 0
    U(i, j) = c(i, j)
   ElseIf i = j Then
    L(i, j) = c(i, j)
    U(i, j) = 1
   ElseIf i > j Then
    L(i, j) = c(i, j)
    U(i, j) = 0
   End If
  Next j
 Next i
End Function
 
'行列の積
Function Prod(X() As Double, Y() As Double, Z() As Double)
 Dim c As Double
 Dim i As Long, j As Long, k As Long, mat1L As Long, mat1U As Long, mat2L As Long, mat2U As Long
 '動的配列をセット
 mat1L = Application.WorksheetFunction.Min(LBound(X, 1), LBound(Y, 1))
 mat1U = Application.WorksheetFunction.Min(UBound(X, 1), UBound(Y, 1))
 mat2L = Application.WorksheetFunction.Min(LBound(X, 2), LBound(Y, 2))
 mat2U = Application.WorksheetFunction.Min(UBound(X, 2), UBound(Y, 2))
 ReDim Z(mat1L To mat1U, mat2L To mat2U)
 '各成分の計算
 For i = 1 To UBound(X, 1)
  For j = 1 To UBound(Y, 2)
  c = 0
   For k = 1 To UBound(X, 2)
    c = c + X(i, k) * Y(k, j)
   Next k
   Z(i, j) = c
  Next j
 Next i
End Function
 
'下三角行列Lの逆行列M
Function LInv(L() As Double, M() As Double)
 Dim c As Double
 Dim i As Long, j As Long, k As Long
 ' 動的配列をセット
 ReDim M(LBound(L, 1) To UBound(L, 1), LBound(L, 2) To UBound(L, 2))
 '各成分の計算
 For i = LBound(L, 1) To UBound(L, 1)
  For j = LBound(L, 2) To UBound(L, 2)
  c = 0
  If i < j Then
   M(i, j) = 0
  ElseIf i = j Then
   M(i, j) = 1 / L(i, j)
  Else
   For k = LBound(L, 2) To UBound(L, 1) - 1
    c = c + L(i, k) * M(k, j)
   Next k
   M(i, j) = -c / L(i, i)
  End If
  Next j
 Next i
End Function

回答
投稿日時: 18/06/06 19:37:32
投稿者: んなっと

小数で計算するとどうしても演算誤差が出ますよね。
そこで、「文字列の分数」に統一するのはどうでしょうか。
+,-,*,/ などを、うしろに追加したTASU(),KAKE(),WARU()に置き換えているだけです。
最大公約数などは、WorksheetFunction.Gcdを使いませんでした。
 
Sub LU分解()
 '記号の形式を定義
  Dim r As Range
  Dim a
  Dim L() As String, U() As String, X() As String, Y() As String, Z() As String, M() As String, N() As String
  Dim i As Long, j As Long
  Set r = Range("C2:F5")
  a = r.Value
  With r.Offset(10, 0).Resize(100)
    .ClearContents
    .NumberFormat = "@"
  End With
 '動的配列をセット
  ReDim L(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  ReDim U(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  'LU分解Call
  Call LU(a, L, U)
  'L行列出力
  Set r = r.Offset(10, 0)
  r.Value = WorksheetFunction.Index(L, 0, 0)
  'U行列出力
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(U, 0, 0)
  'LUが元に戻るかチェック
  Call Prod(L, U, Z)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(Z, 0, 0)
  '下三角行列Lの逆行列M
  Call LInv(L, M)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(M, 0, 0)
  'LMの積が単位行列に戻るかどうかチェック
  Call Prod(L, M, Z)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(Z, 0, 0)
End Sub
'LU分解
Function LU(a, L() As String, U() As String)
  Dim c() As String
  Dim i As Long, j As Long, k As Long
  'L,Uに分配前の成分を計算(※成分が0,1,aijのものは後で配分する)
  ' 動的配列をセット
  ReDim c(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  For i = 1 To UBound(a, 1)
    For j = 1 To UBound(a, 2)
      c(i, j) = a(i, j)
      If j > 1 Then
        If i >= j Then
          For k = 1 To j - 1
            c(i, j) = Tasu(c(i, j), Kake("-1", Kake(c(i, k), c(k, j))))
          Next k
        Else
          For k = 1 To i - 1
            c(i, j) = Tasu(c(i, j), Kake("-1", Kake(c(i, k), c(k, j))))
          Next k
          c(i, j) = Waru(c(i, j), c(i, i))
        End If
      End If
   Next j
  Next i
  'LU行列に分配
  For i = 1 To UBound(a, 1)
    For j = 1 To UBound(a, 2)
      If i < j Then
         L(i, j) = 0
         U(i, j) = c(i, j)
      ElseIf i = j Then
         L(i, j) = c(i, j)
         U(i, j) = 1
      ElseIf i > j Then
         L(i, j) = c(i, j)
         U(i, j) = 0
      End If
    Next j
  Next i
End Function
'行列の積Z
Function Prod(X() As String, Y() As String, Z() As String)
  Dim c As String
  Dim i As Long, j As Long, k As Long, mat1L As Long, mat1U As Long, mat2L As Long, mat2U As Long
  '動的配列をセット
  mat1L = Application.WorksheetFunction.Min(LBound(X, 1), LBound(Y, 1))
  mat1U = Application.WorksheetFunction.Min(UBound(X, 1), UBound(Y, 1))
  mat2L = Application.WorksheetFunction.Min(LBound(X, 2), LBound(Y, 2))
  mat2U = Application.WorksheetFunction.Min(UBound(X, 2), UBound(Y, 2))
  ReDim Z(mat1L To mat1U, mat2L To mat2U)
  '各成分の計算
  For i = 1 To UBound(X, 1)
    For j = 1 To UBound(Y, 2)
      c = "0"
      For k = 1 To UBound(X, 2)
        c = Tasu(c, Kake(X(i, k), Y(k, j)))
      Next k
      Z(i, j) = c
    Next j
  Next i
End Function
'下三角行列Lの逆行列M
Function LInv(L() As String, M() As String)
  Dim c As String
  Dim i As Long, j As Long, k As Long
  ' 動的配列をセット
  ReDim M(LBound(L, 1) To UBound(L, 1), LBound(L, 2) To UBound(L, 2))
  '各成分の計算
  For i = LBound(L, 1) To UBound(L, 1)
    For j = LBound(L, 2) To UBound(L, 2)
      c = "0"
      If i < j Then
        M(i, j) = 0
      ElseIf i = j Then
        M(i, j) = Waru("1", L(i, j))
      Else
        For k = LBound(L, 2) To UBound(L, 1) - 1
          c = Tasu(c, Kake(L(i, k), M(k, j)))
        Next k
        M(i, j) = Kake("-1", Waru(c, L(i, i)))
      End If
    Next j
  Next i
End Function
Function Tasu(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim Var1, Var2
  Dim Bunnbo As String, Bunnsi As String
  If Not Str1 Like "*/*" Then Str1 = Str1 & "/1"
  If Not Str2 Like "*/*" Then Str2 = Str2 & "/1"
  Var1 = Split(Str1, "/")
  Var2 = Split(Str2, "/")
  Bunnbo = Koubaisuu(Var1(1), Var2(1))
  Bunnsi = Var1(0) * Bunnbo / Var1(1) + Var2(0) * Bunnbo / Var2(1)
  Tasu = Yakubunn(Bunnsi, Bunnbo)
End Function
Function Kake(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim Var1, Var2
  If Len(Str1) = 0 Or Len(Str2) = 0 Then
    Kake = "0"
    Exit Function
  End If
  If Not Str1 Like "*/*" Then Str1 = Str1 & "/1"
  If Not Str2 Like "*/*" Then Str2 = Str2 & "/1"
  Var1 = Split(Str1, "/")
  Var2 = Split(Str2, "/")
  Kake = Yakubunn(Var1(0) * Var2(0), Var1(1) * Var2(1))
End Function
Function Waru(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim tmpStr As String
  Dim Var2
  tmpStr = Str2
  If Not tmpStr Like "*/*" Then tmpStr = tmpStr & "/1"
  Var2 = Split(tmpStr, "/")
  tmpStr = Var2(1) & "/" & Var2(0)
  Waru = Kake(Str1, tmpStr)
End Function
Function Yakubunn(ByVal X As String, ByVal Y As String) As String
  If Len(X) = 0 Then Exit Function
  Dim Num As String
  Num = Kouyaku(X, Y)
  Yakubunn = IIf(X * Y < 0, "-", "") & (Abs(X) \ Num)
  If Val(Abs(Y) \ Num) > 1 Then
    Yakubunn = Yakubunn & "/" & (Abs(Y) \ Num)
  End If
End Function
Function Kouyaku(ByVal X As String, ByVal Y As String) As String
  If X Mod Y = 0 Then
    Kouyaku = Abs(Y)
  Else
    Kouyaku = Kouyaku(Y, X Mod Y)
  End If
End Function
Function Koubaisuu(ByVal X As String, ByVal Y As String) As String
  Koubaisuu = X * Y / Kouyaku(X, Y)
End Function
 
5×5行列の時は
  Set r = Range("C2:G6")
に変更。

回答
投稿日時: 18/06/06 21:06:36
投稿者: んなっと

最初のコードのように小数のままで処理して、「見た目だけ」ほぼ正確な分数を表示させたいときは
 
結果を書き込む範囲の表示形式を ユーザー定義
?/????
にしましょう。
 
ただし、整数もすべて 0/1 や 6/1 の形で表示されてしまいます。

回答
投稿日時: 18/06/07 12:34:17
投稿者: んなっと

オーバーフローを減らすために
  Bunnsi = Var1(0) * Bunnbo / Var1(1) + Var2(0) * Bunnbo / Var2(1)
          ↓
  Bunnsi = Var1(0) * (Bunnbo / Var1(1)) + Var2(0) * (Bunnbo / Var2(1))
に修正。
 
途中で行列cの対角成分が0になるの時のエラー対策、「絶対値を比較しての入れ替え」は、ご自分で追加してください。

投稿日時: 18/06/07 12:41:49
投稿者: daipon

お返事ありがとうございます。
少数は2進法への変換によって誤差がでるのですね。
プログラムのところは提案頂いた方法を参考に勉強してみます。
ちなみに2点教えてもらいたいのですが、
 
@最終的に文字列を数値に変換する方法があるでしょうか。
Avbaではなくエクセルで計算すると、少数の誤差がなくゼロになるようですが、内部のプログラムに工夫があるのでしょうか?それとも桁数に制限があってゼロに見えているだけでしょうか。

回答
投稿日時: 18/06/07 13:50:07
投稿者: んなっと

何段階かの小数計算を繰り返せば、演算誤差が拡大します。
あとは、それをどの程度許容するかの問題だけですね。
 
L,U,Mなどはあくまで「文字列の分数」で成分を格納し、
それらをワークシート上に書き込むときだけ小数の「近似値」にする方法。
 
Sub Kinji(Rng As Range)
  Dim rr As Range
  For Each rr In Rng
    rr.Value = Evaluate(rr.Value)
  Next
End Sub
Sub LU分解()
 '記号の形式を定義
  Dim r As Range
  Dim a
  Dim L() As String, U() As String, X() As String, Y() As String, Z() As String, M() As String, N() As String
  Dim i As Long, j As Long
  Set r = Range("C2:F5")
  a = r.Value
  With r.Offset(10, 0).Resize(100)
    .ClearContents
    .NumberFormat = "@"
  End With
 '動的配列をセット
  ReDim L(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  ReDim U(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  'LU分解Call
  Call LU(a, L, U)
  'L行列出力
  Set r = r.Offset(10, 0)
  r.Value = WorksheetFunction.Index(L, 0, 0)
  Kinji r
  'U行列出力
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(U, 0, 0)
  Kinji r
  'LUが元に戻るかチェック
  Call Prod(L, U, Z)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(Z, 0, 0)
  Kinji r
  '下三角行列Lの逆行列M
  Call LInv(L, M)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(M, 0, 0)
  Kinji r
  'LMの積が単位行列に戻るかどうかチェック
  Call Prod(L, M, Z)
  Set r = r.Offset(r.Rows.Count + 1, 0)
  r.Value = WorksheetFunction.Index(Z, 0, 0)
  Kinji r
End Sub
'LU分解
Function LU(a, L() As String, U() As String)
  Dim c() As String
  Dim i As Long, j As Long, k As Long
  'L,Uに分配前の成分を計算(※成分が0,1,aijのものは後で配分する)
  ' 動的配列をセット
  ReDim c(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
  For i = 1 To UBound(a, 1)
    For j = 1 To UBound(a, 2)
      c(i, j) = a(i, j)
      If j > 1 Then
        If i >= j Then
          For k = 1 To j - 1
            c(i, j) = Tasu(c(i, j), Kake("-1", Kake(c(i, k), c(k, j))))
          Next k
        Else
          For k = 1 To i - 1
            c(i, j) = Tasu(c(i, j), Kake("-1", Kake(c(i, k), c(k, j))))
          Next k
          c(i, j) = Waru(c(i, j), c(i, i))
        End If
      End If
   Next j
  Next i
  'LU行列に分配
  For i = 1 To UBound(a, 1)
    For j = 1 To UBound(a, 2)
      If i < j Then
         L(i, j) = 0
         U(i, j) = c(i, j)
      ElseIf i = j Then
         L(i, j) = c(i, j)
         U(i, j) = 1
      ElseIf i > j Then
         L(i, j) = c(i, j)
         U(i, j) = 0
      End If
    Next j
  Next i
End Function
'行列の積Z
Function Prod(X() As String, Y() As String, Z() As String)
  Dim c As String
  Dim i As Long, j As Long, k As Long, mat1L As Long, mat1U As Long, mat2L As Long, mat2U As Long
  '動的配列をセット
  mat1L = Application.WorksheetFunction.Min(LBound(X, 1), LBound(Y, 1))
  mat1U = Application.WorksheetFunction.Min(UBound(X, 1), UBound(Y, 1))
  mat2L = Application.WorksheetFunction.Min(LBound(X, 2), LBound(Y, 2))
  mat2U = Application.WorksheetFunction.Min(UBound(X, 2), UBound(Y, 2))
  ReDim Z(mat1L To mat1U, mat2L To mat2U)
  '各成分の計算
  For i = 1 To UBound(X, 1)
    For j = 1 To UBound(Y, 2)
      c = "0"
      For k = 1 To UBound(X, 2)
        c = Tasu(c, Kake(X(i, k), Y(k, j)))
      Next k
      Z(i, j) = c
    Next j
  Next i
End Function
'下三角行列Lの逆行列M
Function LInv(L() As String, M() As String)
  Dim c As String
  Dim i As Long, j As Long, k As Long
  ' 動的配列をセット
  ReDim M(LBound(L, 1) To UBound(L, 1), LBound(L, 2) To UBound(L, 2))
  '各成分の計算
  For i = LBound(L, 1) To UBound(L, 1)
    For j = LBound(L, 2) To UBound(L, 2)
      c = "0"
      If i < j Then
        M(i, j) = 0
      ElseIf i = j Then
        M(i, j) = Waru("1", L(i, j))
      Else
        For k = LBound(L, 2) To UBound(L, 1) - 1
          c = Tasu(c, Kake(L(i, k), M(k, j)))
        Next k
        M(i, j) = Kake("-1", Waru(c, L(i, i)))
      End If
    Next j
  Next i
End Function
Function Tasu(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim Var1, Var2
  Dim Bunnbo As String, Bunnsi As String
  If Not Str1 Like "*/*" Then Str1 = Str1 & "/1"
  If Not Str2 Like "*/*" Then Str2 = Str2 & "/1"
  Var1 = Split(Str1, "/")
  Var2 = Split(Str2, "/")
  Bunnbo = Koubaisuu(Var1(1), Var2(1))
  Bunnsi = Var1(0) * (Bunnbo / Var1(1)) + Var2(0) * (Bunnbo / Var2(1))
  Tasu = Yakubunn(Bunnsi, Bunnbo)
End Function
Function Kake(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim Var1, Var2
  If Len(Str1) = 0 Or Len(Str2) = 0 Then
    Kake = "0"
    Exit Function
  End If
  If Not Str1 Like "*/*" Then Str1 = Str1 & "/1"
  If Not Str2 Like "*/*" Then Str2 = Str2 & "/1"
  Var1 = Split(Str1, "/")
  Var2 = Split(Str2, "/")
  Kake = Yakubunn(Var1(0) * Var2(0), Var1(1) * Var2(1))
End Function
Function Waru(ByVal Str1 As String, ByVal Str2 As String) As String
  Dim tmpStr As String
  Dim Var2
  tmpStr = Str2
  If Not tmpStr Like "*/*" Then tmpStr = tmpStr & "/1"
  Var2 = Split(tmpStr, "/")
  tmpStr = Var2(1) & "/" & Var2(0)
  Waru = Kake(Str1, tmpStr)
End Function
Function Yakubunn(ByVal X As String, ByVal Y As String) As String
  If Len(X) = 0 Then Exit Function
  Dim Num As String
  Num = Kouyaku(X, Y)
  Yakubunn = IIf(X * Y < 0, "-", "") & (Abs(X) \ Num)
  If Val(Abs(Y) \ Num) > 1 Then
    Yakubunn = Yakubunn & "/" & (Abs(Y) \ Num)
  End If
End Function
Function Kouyaku(ByVal X As String, ByVal Y As String) As String
  X = Abs(X)
  Y = Abs(Y)
  If X Mod Y = 0 Then
    Kouyaku = Y
  Else
    Kouyaku = Kouyaku(Y, X Mod Y)
  End If
End Function
Function Koubaisuu(ByVal X As String, ByVal Y As String) As String
  Koubaisuu = X * Y / Kouyaku(X, Y)
End Function

投稿日時: 18/06/08 12:42:01
投稿者: daipon

お返事ありがとうございます。
参考に勉強してみます。
またお願いします。