lemodule.bas

求解非线性方程的方法

        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

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MatrixToString
'  功能：  将矩阵转换为显示字符串
'  参数：  m   - Integer型变量，矩阵的行数
'          n   - Integer型变量，矩阵的列数
'          mtxA  - Double型m x n二维数组，存放相加的左边矩阵
'          sFormat - 显示矩阵各元素的格式控制字符串
'  返回值：String型，显示矩阵的字符串
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MatrixToString(m As Integer, n As Integer, mtxA() As Double, sFormat As String) As String
    Dim i As Integer, j As Integer
    Dim s As String

    s = ""
    For i = 1 To m
        For j = 1 To n
            s = s + Format(mtxA(i, j), sFormat) + "  "
        Next j
        s = s + Chr(13)
    Next i
    
    MatrixToString = s

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MMul
'  功能：  计算矩阵的乘法
'  参数：  m   - Integer型变量，相乘的左边矩阵的行数
'          n   - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数
'          l   -  Integer型变量，相乘的右边矩阵的列数
'          mtxA  - Double型m x n二维数组，存放相乘的左边矩阵
'          mtxB  - Double型n x l二维数组，存放相乘的右边矩阵
'          mtxC  - Double型m x l二维数组，返回矩阵乘积矩阵
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub MMul(m As Integer, n As Integer, l As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double)
    Dim i As Integer, j As Integer, k As Integer

    For i = 1 To m
        For j = 1 To l
            mtxC(i, j) = 0#
            For k = 1 To n
                mtxC(i, j) = mtxC(i, j) + mtxA(i, k) * mtxB(k, j)
            Next k
        Next j
    Next i

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

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MMqr
'  功能：  用豪斯荷尔德变换法对矩阵进行QR分解
'  参数：   m    - Integer型变量。矩阵的行数， m>=n
'           n    - Integer型变量。矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积维n x n。存放待分解矩阵；返回时，存放分解式中的R矩阵.
'          dblQ  - Double型二维数组，体积为m x m。返回时，存放分解式中的Q矩阵
'  返回值： Boolean型。False，失败；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MMqr(m As Integer, n As Integer, dblA() As Double, dblQ() As Double) As Boolean
    Dim i As Integer, j As Integer, k As Integer, nn As Integer, jj As Integer
    Dim u As Double, alpha As Double, w As Double, t As Double

    If (m < n) Then
        MMqr = False
        Exit Function
    End If

    For i = 1 To m
        For j = 1 To m
            dblQ(i, j) = 0#
            If (i = j) Then
                dblQ(i, j) = 1#
            End If
        Next j
    Next i

    nn = n
    If (m = n) Then
        nn = m - 1
    End If

    For k = 1 To nn
        u = 0#
        For i = k To m
            w = Abs(dblA(i, k))
            If (w > u) Then
                u = w
            End If
        Next i

        alpha = 0#
        For i = k To m
            t = dblA(i, k) / u
            alpha = alpha + t * t
        Next i

        If (dblA(k, k) > 0#) Then
            u = -u
        End If

        alpha = u * Sqr(alpha)
        If (Abs(alpha) + 1# = 1#) Then
            MMqr = False
            Exit Function
        End If

        u = Sqr(2# * alpha * (alpha - dblA(k, k)))
        If ((u + 1#) <> 1#) Then
            dblA(k, k) = (dblA(k, k) - alpha) / u
            For i = k + 1 To m
                dblA(i, k) = dblA(i, k) / u
            Next i

            For j = 1 To m
                t = 0#
                For jj = k To m
                  t = t + dblA(jj, k) * dblQ(jj, j)
                Next jj
                For i = k To m
                    dblQ(i, j) = dblQ(i, j) - 2# * t * dblA(i, k)
                Next i
            Next j

            For j = k + 1 To n
                t = 0#
                For jj = k To m
                  t = t + dblA(jj, k) * dblA(jj, j)
                Next jj
                For i = k To m
                    dblA(i, j) = dblA(i, j) - 2# * t * dblA(i, k)
                Next i
            Next j

            dblA(k, k) = alpha
            For i = k + 1 To m
              dblA(i, k) = 0#
            Next i
        End If
    Next k

    For i = 1 To m - 1
        For j = i + 1 To m
            t = dblQ(i, j)
            dblQ(i, j) = dblQ(j, i)
            dblQ(j, i) = t
        Next j
    Next i

    MMqr = 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

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MUav
'  功能：  用豪斯荷尔德变换及变形QR算法对矩阵进行奇异值分解
'  参数：   m    - Integer型变量。系数矩阵的行数， m>=n
'           n    - Integer型变量。系数矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积维m x n。存放待分解矩阵；
'                  返回时，其对角线存放矩阵的奇异值(以非递增次序排列)，其余元素为0。
'          dblU  - Double型二维数组，体积维m x m。返回时，存放奇异值分解式中的左奇异向量U。
'          dblV  - Double型二维数组，体积维n x n。返回时，存放奇异值分解式中的右奇异向量VT。
'           ka  - Integer型变量。ka=max(m,n)+1
'          eps  - Double型变量。奇异值分解函数中的控制精度参数。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MUav(m As Integer, n As Integer, dblA() As Double, dblU() As Double, dblV() As Double, ka As Integer, eps As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, l As Integer, it As Integer
    Dim ll As Integer, kk As Integer, mm As Integer, nn As Integer, m1 As Integer, ks As Integer
    Dim d As Double, dd As Double, t As Double, sm As Double, sm1 As Double, em1 As Double, sk As Double, ek As Double
    Dim b As Double, c As Double, shh As Double, fg(2) As Double, cs(2) As Double
    ReDim s(ka) As Double, e(ka) As Double, w(ka) As Double

    it = 60
    k = n
    
    If (m - 1 < n) Then
        k = m - 1
    End If

    l = m
    If (n - 2 < m) Then
        l = n - 2
    End If
    If (l < 0) Then
        l = 0
    End If

    ll = k
    
    If (l > k) Then
        ll = l
    End If

    If (ll >= 1) Then
        For kk = 1 To ll
            If (kk <= k) Then
                d = 0#
                For i = kk To m
                    d = d + dblA(i, kk) * dblA(i, kk)
                Next i

                s(kk) = Sqr(d)
                If s(kk) <> 0# Then
                    If (dblA(kk, kk) <> 0#) Then
                        s(kk) = Abs(s(kk))
                        If (dblA(kk, kk) < 0#) Then
                            s(kk) = -s(kk)
                        End If
                    End If
                    For i = kk To m
                        dblA(i, kk) = dblA(i, kk) / s(kk)
                    Next i
                    dblA(kk, kk) = 1# + dblA(kk, kk)
                End If
                s(kk) = -s(kk)
            End If

            If (n >= kk + 1) Then
                For j = kk + 1 To n
                    If ((kk <= k) And (s(kk) <> 0#)) Then
                        d = 0#
                        For i = kk To m
                            d = d + dblA(i, kk) * dblA(i, j)
                        Next i
                        d = -d / dblA(kk, kk)
                        For i = kk To m
                            dblA(i, j) = dblA(i, j) + d * dblA(i, kk)
                        Next i
                    End If
                    e(j) = dblA(kk, j)
                Next j
            End If

            If (kk <= k) Then
                For i = kk To m
                    dblU(i, kk) = dblA(i, kk)
                Next i
            End If