nlmodule.bas

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

        g1 = u * u + v * v
        If (g1 >= g) Then
            If (nIs <> 0) Then
                it = 1
                Call g65c(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
                If (it = 0) Then GoTo l40
            Else
                Call g60c(t, x, y, x1, y1, dx, dy, p, q, k, it)
                If (t >= 0.001) Then GoTo l40
                If (g > 1E-18) Then
                    it = 0
                    Call g65c(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
                    If (it = 0) Then GoTo l40
                End If
            End If

            Call g90c(xr, xi, ar, ai, x, y, p, w, k)
        Else
            g = g1
            x = x1
            y = y1
            nIs = 0
            If (g <= 1E-22) Then
                Call g90c(xr, xi, ar, ai, x, y, p, w, k)
            Else
                u1 = k * ar(1)
                v1 = ai(1)
                For i = 3 To k
                    p = u1 * x
                    q = v1 * y
                    pq = (u1 + v1) * (x + y)
                    u1 = p - q + (k - i + 1) * ar(i - 1)
                    v1 = pq - p - q + (k - i + 1) * ai(i - 1)
                Next i

                p = u1 * u1 + v1 * v1
                If (p <= 1E-20) Then
                    it = 0
                    Call g65c(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
                    If (it = 0) Then GoTo l40
                    Call g90c(xr, xi, ar, ai, x, y, p, w, k)
                Else
                    dx = (u * u1 + v * v1) / p
                    dy = (u1 * v - v1 * u) / p
                    t = 1# + 4# / k
                    Call g60c(t, x, y, x1, y1, dx, dy, p, q, k, it)
                    If (t >= 0.001) Then GoTo l40
                    If (g > 1E-18) Then
                        it = 0
                        Call g65c(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
                        If (it = 0) Then GoTo l40
                    End If
                    Call g90c(xr, xi, ar, ai, x, y, p, w, k)
                End If
            End If
        End If

        If (k = 2) Then
            jt = 0
        Else
            jt = 1
        End If
    Wend

    ' 迭代结束
    NLNdhcRoot = True
   
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：g60c
'  功能：  内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g60c(t As Double, x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double, p As Double, q As Double, k As Integer, it As Integer)
    
    it = 1
    While (it = 1)
        t = t / 1.67
        it = 0
        x1 = x - t * dx
        y1 = y - t * dy
        If (k >= 30) Then
            p = Sqr(x1 * x1 + y1 * y1)
            q = Exp(75# / k)
            If (p >= q) Then it = 1
        End If
    Wend

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：g90c
'  功能：  内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g90c(xr() As Double, xi() As Double, ar() As Double, ai() As Double, x As Double, y As Double, p As Double, w As Double, k As Integer)
    
    Dim i As Integer

    For i = 2 To k
        ar(i) = ar(i) + ar(i - 1) * x - ai(i - 1) * y
        ai(i) = ai(i) + ar(i - 1) * y + ai(i - 1) * x
    Next i

    xr(k - 1) = x * w
    xi(k - 1) = y * w
    k = k - 1
    If (k = 2) Then
        p = ar(1) * ar(1) + ai(1) * ai(1)
        xr(1) = -w * (ar(1) * ar(2) + ai(1) * ai(2)) / p
        xi(1) = w * (ar(2) * ai(1) - ar(1) * ai(2)) / p
    End If
    
End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：g65c
'  功能：  内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g65c(x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double, dd As Double, dc As Double, c As Double, k As Integer, nIs As Integer, it As Integer)
  
    If (it = 0) Then
        nIs = 1
        dd = Sqr(dx * dx + dy * dy)
        If (dd > 1#) Then dd = 1#
        dc = 6.28 / (4.5 * k)
        c = 0#
    End If

    While (True)
        c = c + dc
        dx = dd * Cos(c)
        dy = dd * Sin(c)
        x1 = x + dx
        y1 = y + dy
        If (c <= 6.29) Then
            it = 0
            Exit Sub
        End If
        dd = dd / 1.67
        If (dd <= 0.0000001) Then
            it = 1
            Exit Sub
        End If
        c = 0#
    Wend

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLGrad
'  功能：  使用梯度法求解非线性方程组的一组解，本函数需要调用计算目标函数值的函数Func，其形式为：
'          Function Func(x() as double) as double
'  参数    n    - Integer型变量，方程的个数，也是未知数的个数
'          x    - Double型一维数组，长度为n，存放一组初值x1, x2, …, xn，返回时存放方程组的一组实根
'         nMaxIt - Integer型变量，最大迭代次数
'         eps   - Double型变量，精度控制参数
'          h    - Double型变量，收敛控制参数，一般取1e-5至1e-10之间
'  返回值：Boolean型，False则求解失败，True则求解成功。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLGrad(n As Integer, x() As Double, nMaxIt As Integer, eps As Double, h As Double) As Boolean
    Dim i As Integer, j As Integer
    Dim f As Double, fh As Double, r As Double, sum As Double
    ReDim hx(n) As Double, df(n) As Double

    For j = 1 To nMaxIt
        For i = 1 To n
            If (x(i) = 0) Then
                hx(i) = h
            Else
                hx(i) = h * x(i)
            End If
        Next i
    
        f = Func(x)
        
        If (f <= eps) Then
            NLGrad = True
            Exit Function
        End If
        
        sum = 0
        
        For i = 1 To n
            x(i) = x(i) + hx(i)
            fh = Func(x)
            df(i) = (fh - f) / hx(i)
            sum = sum + df(i) * df(i)
            x(i) = x(i) - hx(i)
        Next i
        
        r = f / sum
    
        For i = 1 To n
            x(i) = x(i) - r * df(i)
        Next i
        
    Next j
    
    NLGrad = False

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLNewtonA
'  功能：  使用拟牛顿法求非线性方程组
'         1）本函数要调用全选主元高斯消去法求解线性代数方程组的函数LEGauss;
'         2）本函数需要调用计算目标函数值的函数Func，其形式为：
'          Sub Func(x() as Double, y() as Double)
'          y(i) 为方程组中各左边多项式与一组x(i)对应的函数值
'  参数    n    - Integer型变量，方程的个数，也是未知数的个数
'          x    - Double型一维数组，长度为n，存放一组初值x1, x2, …, xn，返回时存放方程组的一组实根
'         nMaxIt - Integer型变量，最大迭代次数
'         eps   - Double型变量，精度控制参数
'          h    - Double型变量，增量初值，在本函数中要被破坏
'          t    - Double型变量，控制h大小的变量，0<t<1
'  返回值：Boolean型，True则求解成功；False则求解失败。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLNewtonA(n As Integer, x() As Double, nMaxIt As Integer, eps As Double, h As Double, t As Double) As Boolean
    Dim i As Integer, j As Integer, l As Integer
    Dim am As Double, z As Double, beta As Double, d As Double
    ReDim y(n) As Double, a(n, n) As Double, b(n) As Double

    ' 最大迭代次数
    l = nMaxIt
    
    '精度控制
    am = 1# + eps
    
    ' 迭代求解
    While (am >= eps)
        
        ' 函数值
        Call Func(x, b)
        
        am = 0#
        For i = 1 To n
            z = Abs(b(i))
            If (z > am) Then am = z
        Next i

        If (am >= eps) Then
            l = l - 1
            
            ' 达到最大迭代次数，精度未达到要求，求解失败，返回
            If (l = 0) Then
                NLNewtonA = False
                Exit Function
            End If

            For j = 1 To n
                z = x(j)
                x(j) = x(j) + h
                
                Call Func(x, y)
                
                For i = 1 To n
                    a(i, j) = y(i)
                Next i
                
                x(j) = z
                
            Next j

            ' 高斯消元失败，求解失败，返回
            If Not LEGauss(n, a, b) Then
                NLNewtonA = False
                Exit Function
            End If

            beta = 1#
            For i = 1 To n
                beta = beta - b(i)
            Next i

            ' 方程求解失败，返回
            If (Abs(beta) + 1# = 1#) Then
                NLNewtonA = False
                Exit Function
            End If

            d = h / beta
            For i = 1 To n
                x(i) = x(i) - d * b(i)
            Next i

            h = t * h
            
       End If
    Wend

    ' 方程求解成功，返回
    NLNewtonA = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLMiv
'  功能：  用最小二乘解的广义逆法求非线性方程组
'          1）本函数要调用求解线性最小乘问题的广义逆法的函数LEMiv；
'          2）本函数需要调用计算目标函数值的函数Func，其形式为：
'             Sub Func(x() as Double, y() as Double)
'             y(i) 为方程组中各左边多项式与一组x(i)对应的函数值
'          3）本函数需要调用计算雅可比矩阵函数，其形式为：
'             Sub FuncMJ(x() as Double, p() as Double)
'             p(i, j) 为方程组中各左边多项式fi的对变元xj的偏导数与一组x(i)对应的函数值
'  参数    m    - Integer型变量，非线性方程组中方程的个数
'          n    - Integer型变量，非线性方程组中未知数的个数
'          x    - Double型一维数组，长度为n，存放一组初值x1, x2, …, xn，要求不全为0，返回时存放方程组的
'                 最小二乘解，当m=n时，即是非线性方程组的解
'         eps1   - Double型变量，最小二乘解的精度控制参数
'         eps2   - Double型变量，奇异值分解的精度控制参数
'         ka     - Integer型变量，ka=max(m, n)+1
'  返回值：Boolean型，True则求解成功；False则求解失败。