lemodule.bas

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

        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

        beta = d / e
        d = 0#
        For k = 1 To n
            d = d + r(k) * r(k)
        Next k

        d = Sqr(d)
        
        ' 求解结束，返回
        If (d < eps) Then
            For k = 1 To n
               dblX(k) = x(k, 1)
            Next k
            Exit Sub
        End If

        For k = 1 To n
           p(k, 1) = r(k) - beta * p(k, 1)
        Next k

        i = i + 1
    Wend

    ' 求解结束，返回
    For k = 1 To n
       dblX(k) = x(k, 1)
    Next k

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEMqr
'  功能：  用豪斯荷尔德变换法求解线性最小二乘问题方程组
'  参数：   m    - Integer型变量。系数矩阵的行数， m>=n
'           n    - Integer型变量。系数矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积维n x n。存放系数矩阵；返回时，存放分解式中的R矩阵.
'          dblB  - Double型一维数组，长度为m。存放方程组右端常数向量；返回时，前n个元素存放方程组的最小二乘解。
'          dblQ  - Double型二维数组，体积为m x m。返回时，存放分解式中的Q矩阵
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function LEMqr(m As Integer, n As Integer, dblA() As Double, dblB() As Double, dblQ() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer
    Dim d As Double
    ReDim c(n) As Double

    ' QR分解失败，返回
    If (Not MMqr(m, n, dblA, dblQ)) Then
        LEMqr = False
        Exit Function
    End If

    For i = 1 To n
        d = 0#
        For j = 1 To m
            d = d + dblQ(j, i) * dblB(j)
        Next j
        c(i) = d
    Next i

    dblB(n) = c(n) / dblA(n, n)
    For i = n - 1 To 1 Step -1
        d = 0#
        For j = i + 1 To n
            d = d + dblA(i, j) * dblB(j)
        Next j
        dblB(i) = (c(i) - d) / dblA(i, i)
    Next i

    ' 求解结束，返回
    LEMqr = True
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEMiv
'  功能：  用广义逆法求解线性最小二乘问题方程组
'  参数：   m    - Integer型变量。系数矩阵的行数， m>=n
'           n    - Integer型变量。系数矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积维n x n。存放超定方程组系数矩阵；
'                  返回时，其对角线存放矩阵的奇异值，其余元素为0。
'          dblB  - Double型一维数组，长度为m。存放超定方程组右端常数向量
'          dblX  - Double型一维数组，长度为n。返回时，存放超定方程组的最小二乘解。
'          dblAP - Double型二维数组，体积维n x m。返回时，存放超定方程组系数矩阵A的广义逆A+。
'          dblU  - Double型二维数组，体积维m x m。返回时，存放超定方程组系数矩阵A的奇异值分解式中的
'                  左奇异向量U。
'          dblV  - Double型二维数组，体积维n x n。返回时，存放超定方程组系数矩阵A的奇异值分解式中的
'                  右奇异向量VT。
'           ka  - Integer型变量。ka=max(m,n)+1
'          eps  - Double型变量。奇异值分解函数中的控制精度参数。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function LEMiv(m As Integer, n As Integer, dblA() As Double, dblB() As Double, dblX() As Double, dblAP() As Double, _
dblU() As Double, dblV() As Double, ka As Integer, eps As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer

    If (Not MInv(m, n, dblA, dblAP, dblU, dblV, ka, eps)) Then
        LEMiv = False
        Exit Function
    End If

    For i = 1 To n
        dblX(i) = 0#
        For j = 1 To m
          dblX(i) = dblX(i) + dblAP(i, j) * dblB(j)
        Next j
    Next i

    LEMiv = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：LEModule.bas
'  函数名：LEMorbid
'  功能：  求解病态方程组
'  参数：   n    - Integer型变量，方程组的阶数。
'          dblA  - Double型二维数组，体积维n x n。存放病态方程组系数矩阵。
'          dblB  - Double型一维数组，长度为n，存放方程组右端常数向量。
'          dblX  - Double型一维数组，长度为n。返回时，存放方程组的解向量。
'          eps  - Double型变量。控制精度参数。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function LEMorbid(n As Integer, dblA() As Double, dblB() As Double, dblX() As Double, eps As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, kk As Integer
    Dim q As Double, qq As Double
    ReDim p(n, n) As Double, r(n) As Double, e(n, 1) As Double, x(n, 1) As Double, xx(n) As Double

    i = 60
    For k = 1 To n
        For j = 1 To n
            p(k, j) = dblA(k, j)
        Next j
    Next k

    For k = 1 To n
        x(k, 1) = dblB(k)
    Next k

    For k = 1 To n
        xx(k) = x(k, 1)
    Next k

    ' 全选主元高斯消去法
    If (Not LEGauss(n, p, xx)) Then
        LEMorbid = False
        Exit Function
    End If

    For k = 1 To n
        x(k, 1) = xx(k)
    Next k
    
    q = 1# + eps
    While (q >= eps)
        If (i = 0) Then
            LEMorbid = False
            Exit Function
        End If

        i = i - 1
        
        ' 矩阵乘法
        Call MMul(n, n, 1, dblA, x, e)
        
        For k = 1 To n
            r(k) = dblB(k) - e(k, 1)
        Next k

        For k = 1 To n
            For j = 1 To n
                p(k, j) = dblA(k, j)
            Next j
        Next k

        ' 全选主元高斯消去法
        If (Not LEGauss(n, p, r)) Then
            LEMorbid = False
            Exit Function
        End If

        q = 0#
        For k = 1 To n
            qq = Abs(r(k)) / (1# + Abs(x(k, 1) + r(k)))
            If (qq > q) Then q = qq
        Next k

        For k = 1 To n
            x(k, 1) = x(k, 1) + r(k)
        Next k
    Wend

    ' 解赋值返回
    For k = 1 To n
        dblX(k) = x(k, 1)
    Next k

    LEMorbid = True

End Function