nlmodule.bas

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

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLMiv(m As Integer, n As Integer, x() As Double, eps1 As Double, eps2 As Double, ka As Integer) As Boolean
    Dim i As Integer, j As Integer, k As Integer, l As Integer, kk As Integer
    Dim jt As Boolean
    Dim y(10) As Double, b(10) As Double, alpha As Double, z As Double, h2 As Double, y1 As Double, y2 As Double, y3 As Double, y0 As Double, h1 As Double
    
    ' 分配内存
    ReDim p(m, n) As Double, d(m) As Double, pp(n, m) As Double, dx(n) As Double, u(m, m) As Double, v(n, n) As Double, w(ka) As Double

    ' 迭代次数设为60
    l = 60
    
    ' 控制初值
    alpha = 1#
    
    While (l > 0)
        ' 调用计算目标函数值的函数
        Call Func(x, d)
        
        ' 调用计算雅可比矩阵函数
        Call FuncMJ(x, p)
        
        ' 调用求解线性最小乘问题的广义逆法的函数LEMiv
        jt = LEMiv(m, n, p, d, dx, pp, u, v, ka, eps2)
        
        ' 求解失败，返回
        If (jt = False) Then
            NLMiv = False
            Exit Function
        End If

        ' 继续迭代计算
        j = 0
        jt = True
        h2 = 0#
        
        While (jt)
            jt = False
            If (j <= 2) Then
                z = alpha + 0.01 * j
            Else
                z = h2
            End If

            For i = 1 To n
                w(i) = x(i) - z * dx(i)
            Next i

            Call Func(w, d)
            
            y1 = 0#
            
            For i = 1 To m
                y1 = y1 + d(i) * d(i)
            Next i

            For i = 1 To n
              w(i) = x(i) - (z + 0.00001) * dx(i)
            Next i

            Call Func(w, d)
            
            y2 = 0#
            For i = 1 To m
                y2 = y2 + d(i) * d(i)
            Next i

            y0 = (y2 - y1) / 0.00001
            
            If (Abs(y0) > 0.0000000001) Then
                h1 = y0
                h2 = z
                If (j = 0) Then
                    y(1) = h1
                    b(1) = h2
                Else
                    y(j + 1) = h1
                    kk = 0
                    k = 0
                    While ((kk = 0) And (k <= j - 1))
                        y3 = h2 - b(k + 1)
                        If (Abs(y3) + 1# = 1#) Then
                            kk = 1
                        Else
                            h2 = (h1 - y(k + 1)) / y3
                        End If

                        k = k + 1
                    Wend

                    b(j + 1) = h2
                    If (kk <> 0) Then b(j + 1) = 1E+35
                    h2 = 0#
                    For k = j - 1 To 0 Step -1
                      h2 = -y(k + 1) / (b(k + 2) + h2)
                    Next k

                    h2 = h2 + b(1)
                End If
                j = j + 1
                If (j <= 7) Then
                    jt = True
                Else
                    z = h2
                End If
           End If
        Wend

        alpha = z
        y1 = 0#
        y2 = 0#
        For i = 1 To n
            dx(i) = -alpha * dx(i)
            x(i) = x(i) + dx(i)
            y1 = y1 + Abs(dx(i))
            y2 = y2 + Abs(x(i))
        Next i

        ' 达到精度要求，求解成功，返回
        If (y1 < eps1 * y2) Then
            NLMiv = True
            Exit Function
        End If

        l = l - 1
    Wend

    ' 迭代60次后仍未达到精度要求，求解失败，返回
    NLMiv = False

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLMtclRoot
'  功能：  用蒙特卡洛求非线性方程一个实根，本函数需要调用计算方程左端函数f(x)值的函数Func，其形式为：
'          Function Func(x As Double) As Double
'  参数    x    - Double型变量，存放初值，返回时存放方程的一个实根
'          b    - Double型变量，均匀分布随机数的端点初值
'          m   - Integer型变量，控制调节b的参数
'         eps   - Double型变量，精度控制参数
'  返回值：无
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub NLMtclRoot(x As Double, b As Double, m As Integer, eps As Double)
    Dim k As Integer
    Dim xx As Double, a As Double, y As Double, x1 As Double, y1 As Double

    ' 初值
    a = b
    k = 1
    xx = x
    
    ' 函数值
    y = Func(xx)
    
    ' 迭代求解
    While (a >= eps)
        
        ' 随机数
        Call Randomize(1)
        x1 = Rnd()
        
        x1 = -a + 2# * a * x1
        x1 = xx + x1
        y1 = Func(x1)
        
        k = k + 1
        
        If (Abs(y1) >= Abs(y)) Then
            If (k > m) Then
                k = 1
                a = a / 2#
            End If
        Else
            k = 1
            xx = x1
            y = y1
            
            ' 精度达到要求，求解结束
            If (Abs(y) < eps) Then
                x = xx
                Exit Sub
            End If
        End If
    Wend

    ' 迭代求解结束
    x = xx
    
End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLMtclcRoot
'  功能：  用蒙特卡洛求非线性方程一个复根，本函数需要调用计算方程左端函数f(x)值的函数Func，其形式为：
'          Function Func(x As Double, y As Double) As Double
'          其返回值为||f(x+jy)||
'  参数    x    - Double型变量，存放初值的实部，返回时存放方程的一个复根的实部
'          y    - Double型变量，存放初值的虚部，返回时存放方程的一个复根的虚部
'          b    - Double型变量，均匀分布随机数的端点初值
'          m   - Integer型变量，控制调节b的参数
'         eps   - Double型变量，精度控制参数
'  返回值：无
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub NLMtclcRoot(x As Double, y As Double, b As Double, m As Integer, eps As Double)
    Dim k As Integer
    Dim xx As Double, yy As Double, a As Double, z As Double
    Dim x1 As Double, y1 As Double, z1 As Double

    ' 初值
    a = b
    k = 1
    xx = x
    yy = y

    ' 函数值
    z = Func(xx, yy)

    ' 迭代求解
    While (a >= eps)
        
        ' 随机数
        Call Randomize(1)
        x1 = -a + 2# * a * Rnd()

        x1 = xx + x1

        ' 随机数
        Call Randomize(1)
        y1 = -a + 2# * a * Rnd()

        y1 = yy + y1
        z1 = Func(x1, y1)
        k = k + 1
        If (z1 >= z) Then
            If (k > m) Then
                k = 1
                a = a / 2#
            End If
        Else
            k = 1
            xx = x1
            yy = y1
            z = z1
            
            ' 精度达到要求，求解结束
            If (z < eps) Then
                x = xx
                y = yy
                Exit Sub
            End If
        End If
    Wend

    ' 迭代求解结束
    x = xx
    y = yy

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：NLModule.bas
'  函数名：NLMtcl
'  功能：  用蒙特卡洛法求非线性方程组一组实根，本函数需要调用计算方程左端函数f(x)值的函数Func，其形式为：
'          Function Func(x As Double) As Double
'          其返回值为Sqr(f1*f1 + f2 * f2 + … + fn*fn)
'  参数    n   - Integer型变量，方程的个数，也是未知数的个数
'          x    - Double型一维数组，长度为n，存放一组初值，返回时存放方程的一组实根
'          b    - Double型变量，均匀分布随机数的端点初值
'          m   - Integer型变量，控制调节b的参数
'         eps   - Double型变量，精度控制参数
'  返回值：无
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub NLMtcl(n As Integer, x() As Double, b As Double, m As Integer, eps As Double)
    Dim k As Integer, i As Integer
    Dim a As Double, z As Double, z1 As Double
    ReDim y(n) As Double

    ' 初值
    a = b
    k = 1
    
    ' 函数值
    z = Func(x)

    ' 迭代求解
    While (a >= eps)
        For i = 1 To n
             ' 随机数
            Call Randomize(1)
            y(i) = -a + 2# * a * Rnd() + x(i)
        Next i

        z1 = Func(y)
        
        k = k + 1
        
        If (z1 >= z) Then
            If (k > m) Then
                k = 1
                a = a / 2#
            End If
        Else
            k = 1
            For i = 1 To n
                x(i) = y(i)
            Next i

            z = z1
            
            ' 精度达到要求，求解结束
            If (z < eps) Then Exit Sub
        
        End If
    Wend
End Sub