module1.bas

来自「三角分解法」· BAS 代码 · 共 92 行

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