lemodule.bas

求解非线性方程的方法

        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->
        d = dblB(k)
        dblB(k) = dblB(nIs)
        dblB(nIs) = d
        'A->
        For i = 1 To n
            d = dblA(i, k)
            dblA(i, k) = dblA(i, nJs(k))
            dblA(i, nJs(k)) = d
        Next i
        
        For j = k + 1 To n
            dblA(k, j) = dblA(k, j) / dblA(k, k)
        Next j

        dblB(k) = dblB(k) / dblA(k, k)
        
        ' 回代
        For j = k + 1 To n
          For i = 1 To n
            If i <> k Then
              dblA(i, j) = dblA(i, j) - dblA(i, k) * dblA(k, j)
            End If
          Next i
        Next j
      
        For i = 1 To n
        If i <> k Then
          dblB(i) = dblB(i) - dblA(i, k) * dblB(k)
        End If
        Next i
    Next k
        
    ' 调整解的次序
    For k = n To 1 Step -1
        d = dblB(k)
        dblB(k) = dblB(nJs(k))
        dblB(nJs(k)) = d
    Next k
    
    ' 求解成功
    LEGgje = True
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LETlvs
'  功能：  用列文逊递推算法求解对称托伯利兹方程组
'  参数    n    - Integer型变量，线性代数方程组的阶数
'         dblT  - Double型一维数组，长度为n ，存放对称T型矩阵中的元素
'         dblB  - Double一维数组，长度为n，存放方程组右端的常数向量
'         dblX  - Double一维数组，长度为n，返回存放方程组的解
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LETlvs(n As Integer, dblT() As Double, dblB() As Double, dblX() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer
    Dim a As Double, beta As Double, q As Double, c As Double, h As Double
    ReDim y(n) As Double, s(n) As Double
    
    a = dblT(1)
    If (Abs(a) + 1# = 1#) Then
        LETlvs = False
        Exit Function
    End If

    y(1) = 1#
    dblX(1) = dblB(1) / a

    For k = 1 To n - 1
        beta = 0#
        q = 0#
        For j = 1 To k
            beta = beta + y(j) * dblT(j + 1)
            q = q + dblX(j) * dblT(k - j + 2)
        Next j

        If (Abs(a) + 1# = 1#) Then
            LETlvs = False
            Exit Function
        End If

        c = -beta / a
        s(1) = c * y(k)
        y(k + 1) = y(k)
        If (k <> 1) Then
            For i = 2 To k
                s(i) = y(i - 1) + c * y(k - i + 1)
            Next i
        End If

        a = a + c * beta
        If (Abs(a) + 1# = 1#) Then
            LETlvs = False
            Exit Function
        End If

        h = (dblB(k + 1) - q) / a
        For i = 1 To k
            dblX(i) = dblX(i) + h * s(i)
            y(i) = s(i)
        Next i

        dblX(k + 1) = h * y(k + 1)
    Next k

    LETlvs = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LESeidel
'  功能：  用高斯-赛德尔迭代求解系数矩阵主对角线占绝对优势线性方程组
'  参数    n    - Integer型变量，线性代数方程组的阶数
'         dblA  - Double型n x n二维数组，存放系数矩阵
'         dblB  - Double一维数组，长度为n，存放方程组右端的常数向量
'         dblX  - Double型一维数组，长度为n。返回方程组的解。
'         eps  - Double型变量。给定的精度要求。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function LESeidel(n As Integer, dblA() As Double, dblB() As Double, dblX() As Double, eps As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer
    Dim p As Double, q As Double, s As Double, t As Double

    ' 校验系数矩阵是否主对角线占绝对优势
    For i = 1 To n
        p = 0#
        dblX(i) = 0#
        For j = 1 To n
            If (i <> j) Then
                p = p + Abs(dblA(i, j))
            End If
        Next j
        If (p >= Abs(dblA(i, i))) Then
            LESeidel = False
            Exit Function
        End If
    Next i

    ' 迭代求解
    p = eps + 1#
    While (p >= eps)
        p = 0#
        For i = 1 To n
            t = dblX(i)
            s = 0#
            For j = 1 To n
                If (j <> i) Then
                    s = s + dblA(i, j) * dblX(j)
                End If
            Next j

            dblX(i) = (dblB(i) - s) / dblA(i, i)
            q = Abs(dblX(i) - t) / (1# + Abs(dblX(i)))
            If (q > p) Then
                p = q
            End If
        Next i
    Wend

    ' 求解成功
    LESeidel = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEGrad
'  功能：  用共轭梯度法是求解n阶对称正定方程组
'  参数    n    - Integer型变量，线性代数方程组的阶数
'         dblA  - Double型n x n二维数组，存放对称正定系数矩阵
'         dblB  - Double一维数组，长度为n，存放方程组右端的常数向量
'         dblX  - Double型一维数组，长度为n。返回方程组的解。
'         eps  - Double型变量。给定的精度要求。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub LEGrad(n As Integer, dblA() As Double, dblB() As Double, dblX() As Double, eps As Double)
    ' 局部变量
    Dim i As Integer, k As Integer
    ReDim p(n, 1) As Double, r(n) As Double, s(n, 1) As Double, q(n, 1) As Double, x(n, 1) As Double
    Dim alpha As Double, beta As Double, d As Double, e As Double

    ' 初始化
    For i = 1 To n
        x(i, 1) = 0#
        p(i, 1) = dblB(i)
        r(i) = dblB(i)
    Next i

    ' 循环求解
    i = 1
    While (i <= n)
    
        ' 矩阵乘法
        Call MMul(n, n, 1, dblA, p, s)
        
        d = 0#
        e = 0#
        
        For k = 1 To n
            d = d + p(k, 1) * dblB(k)
            e = e + p(k, 1) * s(k, 1)
        Next k

        alpha = d / e
        For k = 1 To n
           x(k, 1) = x(k, 1) + alpha * p(k, 1)
        Next k

        ' 矩阵乘法
        Call MMul(n, n, 1, dblA, x, q)

        d = 0#
        For k = 1 To n
            r(k) = dblB(k) - q(k, 1)
            d = d + r(k) * s(k, 1)
        Next k