lemodule.bas

科学与工程数值计算算法（Visual Basic版） 附赠的光盘包含了本书中全部的源代码

            q = -dblAI(k, j) * dblAI(k, k)
            s = (dblAR(k, k) - dblAI(k, k)) * (dblAR(k, j) + dblAI(k, j))
            dblAR(k, j) = (p - q) / d
            dblAI(k, j) = (s - p - q) / d
        Next j
        For j = 1 To m
            p = dblBR(k, j) * dblAR(k, k)
            q = -dblBI(k, j) * dblAI(k, k)
            s = (dblAR(k, k) - dblAI(k, k)) * (dblBR(k, j) + dblBI(k, j))
            dblBR(k, j) = (p - q) / d
            dblBI(k, j) = (s - p - q) / d
        Next j
        
        For i = 1 To n
            If i <> k Then
                For j = k + 1 To n
                    p = dblAR(i, k) * dblAR(k, j)
                    q = dblAI(i, k) * dblAI(k, j)
                    s = (dblAR(i, k) + dblAI(i, k)) * (dblAR(k, j) + dblAI(k, j))
                    dblAR(i, j) = dblAR(i, j) - p + q
                    dblAI(i, j) = dblAI(i, j) - s + p + q
                Next j
                For j = 1 To m
                    p = dblAR(i, k) * dblBR(k, j)
                    q = dblAI(i, k) * dblBI(k, j)
                    s = (dblAR(i, k) + dblAI(i, k)) * (dblBR(k, j) + dblBI(k, j))
                    dblBR(i, j) = dblBR(i, j) - p + q
                    dblBI(i, j) = dblBI(i, j) - s + p + q
                Next j
            End If
        Next i
    Next k
        
    ' 调整解的次序
    For k = n To 1 Step -1
        If nJs(k) <> k Then
            For j = 1 To m
                p = dblBR(k, j)
                dblBR(k, j) = dblBR(nJs(k), j)
                dblBR(nJs(k), j) = p
                
                p = dblBI(k, j)
                dblBI(k, j) = dblBI(nJs(k), j)
                dblBI(nJs(k), j) = p
            Next j
        End If
    Next k
    
    ' 求解成功
    LECpxGaussJordan = True
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LETrid
'  功能：  使用追赶法求解三对角线线性代数方程组
'  参数    n   - Integer型变量，线性代数方程组的阶数
'          m   - Integer型变量，n阶三对角线矩阵三对角线上元素的个数，即数组b的长度。
'                 它的值应为m = 3n -2。函数应对此值进行检验。
'         dblB － Double型一维数组，长度为m。以行为主存放三对角线矩阵中三对角线上的元素，即b中依次存放下列元素：
'                  a11，a12，a21，a22，a23，a32，a33，a34，…，an,n-1，an,n
'         dblD － Double型一维数组，长度为n。作为传入参数，存放方程组右端的常数向量。
'                  函数返回时，此数组中存放着方程组的解向量。
'  返回值： Integer型。小于0，m的值不正确；为0，求解失败，无解；大于0，成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LETrid(n As Integer, m As Integer, dblB() As Double, dblD() As Double) As Integer
    ' 局部变量
    Dim k As Integer, j As Integer
    Dim s As Double

    ' 参数校验
    If (m <> (3 * n - 2)) Then
        LETrid = -1
        Exit Function
    End If

    ' 求解
    For k = 1 To n - 1
        j = 3 * (k - 1) + 1
        s = dblB(j)
        
        ' 无解，返回
        If (Abs(s) + 1# = 1#) Then
            LETrid = 0
            Exit Function
        End If

        dblB(j + 1) = dblB(j + 1) / s
        dblD(k) = dblD(k) / s
        dblB(j + 3) = dblB(j + 3) - dblB(j + 2) * dblB(j + 1)
        dblD(k + 1) = dblD(k + 1) - dblB(j + 2) * dblD(k)
    Next k

    s = dblB(3 * n - 2)
     
     ' 无解，返回
    If (Abs(s) + 1# = 1#) Then
        LETrid = 0
        Exit Function
    End If

    dblD(n) = dblD(n) / s
    For k = n - 1 To 1 Step -1
      dblD(k) = dblD(k) - dblB(3 * (k - 1) + 2) * dblD(k + 1)
    Next k

    ' 求解成功
    LETrid = 1

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEBand
'  功能：  一般带型方程组的求解
'  参数    n   - Integer型变量，线性代数方程组的阶数
'          m  - Integer型变量，为方程组右端的常数向量的个数
'          l   - Integer型变量，为系数矩阵的半带宽。
'          il   - Integer型变量，为系数矩阵的带宽。
'         dblB - Double型n x il二维数组，存放带型矩阵A中带区内的元素
'         dblD - Double型n x m二维数组，作为传入参数，存放方程组右端的m组常数向量。
'                函数返回时，其中存放着m组解向量。
'  返回值： Integer型。小于0，参数中半带宽l与带宽il的关系不对；为0，系数矩阵A奇异，无解；大于0成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LEBand(n As Integer, m As Integer, l As Integer, il As Integer, dblB() As Double, dblD() As Double) As Integer
    ' 局部变量
    Dim ls As Integer, k As Integer, i As Integer, j As Integer, nIs As Integer
    Dim p As Double, t As Double

    ' 参数校验
    If (il <> (2 * l + 1)) Then
        LEBand = -1
        Exit Function
    End If

    ' 求解
    ls = l
    For k = 1 To n - 1
        p = 0#
        For i = k To ls
            t = Abs(dblB(i, 1))
            If t > p Then
                p = t
                nIs = i
            End If
        Next i

        ' 无解，返回
        If (p + 1# = 1#) Then
            LEBand = 0
            Exit Function
        End If

        For j = 1 To m
            t = dblD(k, j)
            dblD(k, j) = dblD(nIs, j)
            dblD(nIs, j) = t
        Next j

        For j = 1 To il
            t = dblB(k, j)
            dblB(k, j) = dblB(nIs, j)
            dblB(nIs, j) = t
        Next j

        For j = 1 To m
          dblD(k, j) = dblD(k, j) / dblB(k, 1)
        Next j

        For j = 2 To il
          dblB(k, j) = dblB(k, j) / dblB(k, 1)
        Next j

        For i = k + 1 To ls + 1
            t = dblB(i, 1)
            For j = 1 To m
                dblD(i, j) = dblD(i, j) - t * dblD(k, j)
            Next j

            For j = 2 To il
                dblB(i, j - 1) = dblB(i, j) - t * dblB(k, j)
            Next j

            dblB(i, il) = 0#
        Next i
        If (ls <> n - 1) Then
            ls = ls + 1
        End If
    Next k

    p = dblB(n, 1)
    
    ' 无解，返回
    If (Abs(p) + 1# = 1#) Then
        LEBand = 0
        Exit Function
    End If

    For j = 1 To m
        dblD(n, j) = dblD(n, j) / p
    Next j

    ls = 1
    For i = n - 1 To 1 Step -1
        For k = 1 To m
            For j = 2 To ls + 1
                dblD(i, k) = dblD(i, k) - dblB(i, j) * dblD(i + j - 1, k)
            Next j
        Next k
        If (ls <> (il - 1)) Then
            ls = ls + 1
        End If
    Next i

    ' 求解成功
    LEBand = 1

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEDjn
'  功能：  用分解法求解对称方程组
'  参数    n   - Integer型变量，线性代数方程组的阶数
'          m  - Integer型变量，为方程组右端的常数向量的个数
'         dblA - Double型n x n二维数组，存放系数矩阵
'         dblC - Double型n x m二维数组，作为传入参数，存放方程组右端的m组常数向量。
'                函数返回时，其中存放着m组解向量。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LEDjn(n As Integer, m As Integer, dblA() As Double, dblC() As Double) As Integer
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, k1 As Integer, k2 As Integer, k3 As Integer
    
    ' 无解，返回
    If (Abs(dblA(1, 1)) + 1# = 1#) Then
        LEDjn = False
        Exit Function
    End If

    For i = 2 To n
      dblA(i, 1) = dblA(i, 1) / dblA(1, 1)
    Next i

    For i = 2 To n - 1
        For j = 2 To i
            dblA(i, i) = dblA(i, i) - dblA(i, j - 1) * dblA(i, j - 1) * dblA(j - 1, j - 1)
        Next j

        ' 无解，返回
        If (Abs(dblA(i, i)) + 1# = 1#) Then
            LEDjn = False
            Exit Function
        End If

        For k = i + 1 To n
            For j = 2 To i
                dblA(k, i) = dblA(k, i) - dblA(k, j - 1) * dblA(i, j - 1) * dblA(j - 1, j - 1)
            Next j
            dblA(k, i) = dblA(k, i) / dblA(i, i)
        Next k
    Next i

    For j = 2 To n
        dblA(n, n) = dblA(n, n) - dblA(n, j - 1) * dblA(n, j - 1) * dblA(j - 1, j - 1)
    Next j
    
    ' 无解，返回
    If (Abs(dblA(n, n)) + 1# = 1#) Then
        LEDjn = False
        Exit Function
    End If

    For j = 1 To m
        For i = 2 To n
            For k = 2 To i
                dblC(i, j) = dblC(i, j) - dblA(i, k - 1) * dblC(k - 1, j)
            Next k
        Next i
    Next j

    For i = 2 To n
        For j = i To n
            dblA(i - 1, j) = dblA(i - 1, i - 1) * dblA(j, i - 1)
        Next j
    Next i

    For j = 1 To m
        dblC(n, j) = dblC(n, j) / dblA(n, n)

        For k = 2 To n
            k1 = n - k + 2
            For k2 = k1 To n
                k3 = n - k + 1
                dblC(k3, j) = dblC(k3, j) - dblA(k3, k2) * dblC(k2, j)
            Next k2
            dblC(k3, j) = dblC(k3, j) / dblA(k3, k3)
        Next k
    Next j

    ' 求解成功
    LEDjn = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LECholesky
'  功能：  用乔里斯基分解法求解正定方程组
'  参数    n   - Integer型变量，线性代数方程组的阶数
'          m  - Integer型变量，为方程组右端的常数向量的个数
'         dblA - Double型n x n二维数组，存放系数矩阵（应为对称正定矩阵）；返回时，其上三角部分存放分解后的U矩阵
'         dblD - Double型n x m二维数组，作为传入参数，存放方程组右端的m组常数向量。
'                函数返回时，其中存放着m组解向量。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LECholesky(n As Integer, m As Integer, dblA() As Double, dblD() As Double) As Integer
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer

    ' 矩阵非正定，求解失败
    If ((dblA(1, 1) + 1# = 1#) Or (dblA(1, 1) < 0#)) Then
        LECholesky = False
        Exit Function
    End If

    dblA(1, 1) = Sqr(dblA(1, 1))
    For j = 2 To n
        dblA(1, j) = dblA(1, j) / dblA(1, 1)
    Next j

    For i = 2 To n
        For j = 2 To i
            dblA(i, i) = dblA(i, i) - dblA(j - 1, i) * dblA(j - 1, i)
        Next j

        ' 求解失败
        If ((dblA(i, i) + 1# = 1#) Or (dblA(i, i) < 0#)) Then
            LECholesky = False
            Exit Function
        End If

        dblA(i, i) = Sqr(dblA(i, i))
        If (i <> n) Then
            For j = i + 1 To n
                For k = 2 To i
                  dblA(i, j) = dblA(i, j) - dblA(k - 1, i) * dblA(k - 1, j)
                Next k
                dblA(i, j) = dblA(i, j) / dblA(i, i)
            Next j
        End If
    Next i

    For j = 1 To m
        dblD(1, j) = dblD(1, j) / dblA(1, 1)

        For i = 2 To n
            For k = 2 To i
                dblD(i, j) = dblD(i, j) - dblA(k - 1, i) * dblD(k - 1, j)
            Next k

            dblD(i, j) = dblD(i, j) / dblA(i, i)
        Next i
    Next j

    For j = 1 To m
        dblD(n, j) = dblD(n, j) / dblA(n, n)

        For k = n To 2 Step -1
            For i = k To n
                dblD(k - 1, j) = dblD(k - 1, j) - dblA(k - 1, i) * dblD(i, j)
            Next i
            dblD(k - 1, j) = dblD(k - 1, j) / dblA(k - 1, k - 1)
        Next k
    Next j

    ' 求解成功
    LECholesky = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEGgje
'  功能：  用全选主元高斯－约当消去法求解稀疏方程组
'  参数    n    - Integer型变量，线性代数方程组的阶数
'         dblA  - Double型n x n二维数组，存放系数矩阵（应为稀疏矩阵）
'         dblB  - Double一维数组，长度为n，存放方程组右端的常数向量；返回时存放方程组的解
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LEGgje(n As Integer, dblA() As Double, dblB() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer
    Dim nIs As Integer
    ReDim nJs(n) As Integer
    Dim d As Double, q As Double
    
    ' 开始求解
    For k = 1 To n
        q = 0#
        
        ' 归一
        For i = k To n
            For j = k To n
                If Abs(dblA(i, j)) > q Then
                    q = Abs(dblA(i, j))
                    nJs(k) = j
                    nIs = i
                End If
            Next j
        Next i
        
        ' 无解，返回
        If q + 1# = 1# Then
            LEGgje = False
            Exit Function
        End If
        
        ' 消元
        ' A->
        For j = k To n
            d = dblA(k, j)
            dblA(k, j) = dblA(nIs, j)
            dblA(nIs, j) = d
        Next j
        ' B->