nlmodule.bas

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

Attribute VB_Name = "NLModule"
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  功能：  求解非线性方程和方程组
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Option Explicit

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLBisectRoot
'  功能：  使用对分法求解非线性方程的实根，本函数需要调用计算方程左端函数f(x)值的函数Func，其形式为：
'          Function Func(x As Double) As Double
'  参数    m    - Integer型变量，在[a, b]内实根个数的预估值
'          a    - Double型变量，求根区间的左端点
'          b    - Double型变量，求根区间的右端点
'          h    - Double型变量，搜索求根时采用的步长
'          x    - Double型一维数组，长度为m。返回在区间[a, b]内搜索到的实根，实根个数由函数值返回
'         eps   - Double型变量，精度控制参数
'  返回值：Integer型，求得的实根的个数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLBisectRoot(m As Integer, a As Double, b As Double, h As Double, x() As Double, eps As Double) As Integer
    Dim n As Integer, js As Integer
    Dim z As Double, y As Double, z1 As Double, y1 As Double, z0 As Double, y0 As Double

    ' 根的个数清0
    n = 0
    
    ' 左边界函数值
    z = a
    y = Func(z)
    
    ' 迭代求解，直到到达右边界
    While ((z <= b + h / 2#) And (n <> m))
        ' 如果精度满足要求，则求得一个实根，继续计算下一步
        If (Abs(y) < eps) Then
            n = n + 1
            x(n) = z
            z = z + h / 2#
            y = Func(z)
        Else
            z1 = z + h
            y1 = Func(z1)
            If (Abs(y1) < eps) Then
                n = n + 1
                x(n) = z1
                z = z1 + h / 2#
                y = Func(z)
            Else
                If (y * y1 > 0#) Then
                    y = y1
                    z = z1
                Else
                    js = 0
                    While (js = 0)
                        If (Abs(z1 - z) < eps) Then
                            n = n + 1
                            x(n) = (z1 + z) / 2#
                            z = z1 + h / 2#
                            y = Func(z)
                            js = 1
                        Else
                            z0 = (z1 + z) / 2#
                            y0 = Func(z0)
                            If (Abs(y0) < eps) Then
                                x(n) = z0
                                n = n + 1
                                js = 1
                                z = z0 + h / 2#
                                y = Func(z)
                            Else
                                If ((y * y0) < 0#) Then
                                    z1 = z0
                                    y1 = y0
                                Else
                                    z = z0
                                    y = y0
                                End If
                            End If
                        End If
                    Wend
                End If
            End If
        End If
    Wend

    ' 返回求得的根的个数
    NLBisectRoot = n

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLNewtonRoot
'  功能：  使用牛顿迭代法求解非线性方程的实根，本函数需要调用计算方程左端函数f(x)值的过程Func，其形式为：
'          Sub Func(x As Double, y() as double)
'          y(1) 返回f(x)的值
'          y(2) 返回f'(x)的值
'  参数    x    - Double型变量，输入时存放迭代初值；返回时存放迭代终值，即方程的根
'         js    - Integer型变量，最大迭代次数
'         eps   - Double型变量，精度控制参数
'  返回值：Integer型，若小于0，则表示在求解失败；若等于最大迭代次数js，则表示迭代了js次还未满足精度要求，
'         返回的实根只作为参考；若大于等于0且小于最大迭代次数js，则表示正常返回。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLNewtonRoot(x As Double, js As Integer, eps As Double) As Integer
    Dim k As Integer, l As Integer
    Dim y(2) As Double, d As Double, p As Double, x0 As Double, x1 As Double

    ' 初值
    l = js
    k = 1
    x0 = x
    
    ' 计算f(x)和f'(x)
    Call Func(x0, y)
    
    d = eps + 1#
    
    ' 迭代计算
    While ((d >= eps) And (l <> 0))
        ' 求解失败，返回
        If (Abs(y(2)) + 1# = 1#) Then
            NLNewtonRoot = -1
            Exit Function
        End If

        x1 = x0 - y(1) / y(2)

        Call Func(x1, y)

        d = Abs(x1 - x0)
        p = Abs(y(1))

        If (p > d) Then d = p
        
        x0 = x1
        l = l - 1
    Wend

    ' 求解结束
    x = x1
    k = js - l

    NLNewtonRoot = k

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLAitkenRoot
'  功能：  使用埃特金迭代法求解非线性方程的一个实根，本函数需要调用计算方程左端函数f(x)值的过程Func，其形式为：
'          Function Func(x as double) as double
'  参数    x    - Double型变量，输入时存放迭代初值；返回时存放迭代终值，即方程的根
'         js    - Integer型变量，最大迭代次数
'         eps   - Double型变量，精度控制参数
'  返回值：Integer型，若为0，则表示迭代了js次还为满足精度要求，返回的实根只作为参考；
'                     若大于0，则表示正常返回。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLAitkenRoot(x As Double, js As Integer, eps As Double) As Integer
    Dim flag As Integer, l As Integer
    Dim u As Double, v As Double, x0 As Double

    ' 迭代初值
    l = 0
    x0 = x
    flag = 0
    
    ' 迭代求解
    While ((flag = 0) And (l <> js))
        l = l + 1
        u = Func(x0)
        v = Func(u)
        If (Abs(u - v) < eps) Then
            x0 = v
            flag = 1
        Else
            x0 = v - (v - u) * (v - u) / (v - 2# * u + x0)
        End If
    Wend

    ' 求解结束
    x = x0
    l = js - l

    NLAitkenRoot = l

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLPqRoot
'  功能：  使用连分式解法求解非线性方程的一个实根，本函数需要调用计算方程左端函数f(x)值的过程Func，其形式为：
'          Function Func(x as double) as double
'  参数    x    - Double型变量，输入时存放迭代初值；返回时存放迭代终值，即方程的根
'         eps   - Double型变量，精度控制参数
'  返回值：Boolean型，若为False，则表示迭代了10次还未满足精度要求，返回的实根只作为参考；
'                     若为True，则表示正常返回。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLPqRoot(x As Double, eps As Double) As Boolean
    Dim i As Integer, j As Integer, m As Integer, it As Integer, l As Integer
    Dim a(10) As Double, y(10) As Double, z As Double, h As Double, x0 As Double, q As Double

    ' 迭代初值
    l = 10
    q = 1E+35
    x0 = x
    h = 0#
    
    ' 迭代计算
    While (l <> 0)
        l = l - 1
        j = 0
        it = l
        
        ' 最多迭代到第7节连分式
        While (j <= 7)
             If (j <= 2) Then
                z = x0 + 0.1 * j
             Else
                z = h
             End If

             ' 函数值
             y(j + 1) = Func(z)
             h = z
             If (j = 0) Then
                a(1) = z
             Else
                 m = 0
                 i = 0
                 While ((m = 0) And (i <= j - 1))
                     If (Abs(h - a(i + 1)) + 1# = 1#) Then
                        m = 1
                     Else
                        h = (y(j + 1) - y(i + 1)) / (h - a(i + 1))
                     End If

                     i = i + 1
                 Wend
                 a(j + 1) = h
                 If (m <> 0) Then a(j + 1) = q
                 h = 0#
                 For i = j To 1 Step -1
                     If (Abs(a(i + 1) + h) + 1# = 1#) Then
                        h = q
                     Else
                        h = -y(i) / (a(i + 1) + h)
                     End If
                 Next i
                 h = h + a(1)
             End If

             ' 精度达到要求，则结束迭代
             If (Abs(y(j + 1)) >= eps) Then
                j = j + 1
             Else
                j = 10
                l = 0
             End If

         Wend
         
         x0 = h
    
    Wend

    ' 求解结束
    x = h
    
    ' 判断解的合理性
    If it = 0 Then
        NLPqRoot = False
        Exit Function
    End If
     
    ' 正常解
    NLPqRoot = True
     
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLQrRoot
'  功能：  使用QR方法求解实系数代数方程的全部根，本函数需要调用用QR方法计算上－H阵全部特征值的函数
'          MhbergEigenv
'  参数    n    - 多项式方程的次数
'          a    - Double型一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
'          xr    - Double型一维数组，长度为n，返回n个根的实部
'          xi    - Double型一维数组，长度为n，返回n个根的虚部
'         eps   - Double型变量，精度控制参数
'         nMaxItNum    - Integer型变量，控制QR方法的最大迭代次数
'  返回值：Boolean型，若为False，则表示在QR方法中迭代已超过最大迭代次数但还未满足精度要求；
'                     若为True，则表示求解成功。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLQrRoot(n As Integer, a() As Double, xr() As Double, xi() As Double, eps As Double, nMaxItNum As Integer) As Boolean
    Dim i As Integer, j As Integer
    ReDim q(n, n) As Double

    ' 用最高幂系数约化其余系数
    For j = 1 To n
      q(1, j) = -a(j + 1) / a(1)
    Next j

    ' 构造上－H阵
    For i = 2 To n
        For j = 1 To n
            q(i, j) = 0#
        Next j
    Next i

    For i = 2 To n
        q(i, i - 1) = 1#
    Next i

    ' 求上－H阵的特征值，即方程的全部根
    NLQrRoot = MHbergEigenv(n, q, xr, xi, eps, nMaxItNum)

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''