newt.bas

来自「NEWT rutine from Numerical Recipes in C 」· BAS 代码 · 共 122 行

Attribute VB_Name = "Module1"
Option Explicit

Const maxits As Integer = 200
Const tolf As Double = 0.0001
Const tolmin As Double = 0.000001
Const tolx As Double = 0.0000001
Const stpmx As Double = 100#

Public nn As Integer
Public fvec() As Double

Public Sub newt(ByRef x() As Double, ByVal n As Integer, ByRef check As Integer)

Dim i As Integer, its As Integer, j As Integer, indx() As Double
Dim d As Double, den As Double, f As Double, fold As Double, stpmax As Double
Dim sum As Double, temp As Double, test As Double, fjac() As Double, g() As Double
Dim p() As Double, xold() As Double

ReDim indx(1 To n)
ReDim fjac(1 To n, 1 To n)
ReDim g(1 To n)
ReDim p(1 To n)
ReDim xold(1 To n)
ReDim fvec(1 To n)

nn = n

nrfuncv = vecfunc
f = fmin(x)
test = 0#
For i = 1 To n
    If Abs(fvec(i) > test) Then
        test = Abs(fvec(i))
    End If
Next i
If test < 0.01 * tolf Then
    check = 0
    Exit Sub
End

sum = 0#
For i = 1 To n
    sum = sum + cuadrado(x(i))
Next i
stpmax = stpmx * fmax(Sqr(sum), CDbl(n))

For its = 1 To maxits
    fdjac n, x, fvec, fjac, vecfunc
    
        For i = 1 To n
            sum = 0#
            For j = 1 To n
                sum = sum + fjac(j, i) * fvec(j)
                g(i) = sum
            Next j
        Next i
        
        For i = 1 To n
            xold(i) = x(i)
        Next i
        fold = f
        
        For i = 1 To n
            p(i) = -fvec(i)
        Next i
        ludcmp fjac, n, indx, d
        lubksb fjac, n, indx, p
        lnsrch n, xold, fold, g, p, x, f, , stpmax, check, fmin
        test = 0#
        For i = 1 To n
            If Abs(fvec(i) > test) Then
                test = Abs(fvec(i))
            End If
        Next i
        If test < tolf Then
            check = 0
            Exit Sub
        End If
        
        If CBool(check) Then
            test = 0#
            den = fmax(f, 0.5 * n)
            For i = 1 To n
                temp = Abs(g(i)) * fmax(Abs(x(i)), 1#) / den
                If temp > test Then
                    test = temp
                End If
            Next i
            If test < tolmin Then
                check = 1
            Else
                check = 0
            End If
            Exit Sub
        End If
        
        test = 0#
        For i = 1 To n
            temp = (Abs(x(i) - xold(i))) / fmax(Abs(x(i)), 1#)
            If temp > test Then
                test = temp
            End If
        Next i
        
        If test < tolx Then
            Exit Sub
        End If
        
Next its
MsgBox ("MAXITS exceeded in newt")
    


End Sub

Public Function cuadrado(ByVal a As Double) As Double

cuadrado = a * a

End Function