module1.bas

三角分解法

Attribute VB_Name = "triangular_decomposition"
Option Explicit

Public Function triangular(n As Integer, a() As Double, b() As Double, x() As Double)
'三角分解法
Dim i, j, k, q, p As Integer
Dim U() As Double
Dim L() As Double
Dim D() As Double
Dim y() As Double
Dim z() As Double
ReDim y(n) As Double
ReDim z(n) As Double
ReDim U(n, n) As Double
ReDim L(n, n) As Double
ReDim D(n, n) As Double
For k = 1 To n
    For j = 1 To k
        If j = k Then
            L(k, j) = 1
        Else
            For p = 1 To j - 1
                L(k, j) = L(k, j) + L(k, p) * U(p, j)
            Next p
            L(k, j) = a(k, j) - L(k, j)
            L(k, j) = L(k, j) / U(j, j)
        End If
    Next j

    For i = 1 To k
        For p = 1 To i - 1
            U(i, k) = U(i, k) + L(i, p) * U(p, k)
        Next p
        U(i, k) = a(i, k) - U(i, k)
    Next i
    
Next k

For i = 1 To n
    For j = i To n
        If i = j Then
            D(i, i) = U(i, j)
            U(i, j) = U(i, j) / D(i, i)
        Else
            U(i, j) = U(i, j) / D(i, i)
        End If
    Next j
Next i
        
For i = 1 To n
    For j = 1 To i - 1
        y(i) = y(i) + y(j) * L(i, j)
    Next j
    y(i) = b(i) - y(i)
Next i

For i = 1 To n
    z(i) = y(i) / D(i, i)
Next i

For i = n To 1 Step -1
    For j = n To i + 1 Step -1
        x(i) = x(i) + x(j) * U(i, j)
    Next j
    x(i) = z(i) - x(i)
Next i

    
       
j = 1
Debug.Print "L="
For i = 1 To n
Debug.Print L(i, j), L(i, j + 1), L(i, j + 2)
Next i
Debug.Print "U="
For i = 1 To n
Debug.Print U(i, j), U(i, j + 1), U(i, j + 2)
Next i
Debug.Print "D="
For i = 1 To n
Debug.Print D(i, j), D(i, j + 1), D(i, j + 2)
Next i

Debug.Print "x="
For i = 1 To n
Debug.Print x(i)
Next i


End Function