d9r8.bas

矩阵特征值的求解过程之一

Attribute VB_Name = "Module1"
    Public NCOM, PCOM(50), XICOM(50)
    Sub LINMIN(P(), XI(), N, FRET)
        TOL = 0.0001
        NCOM = N
        For J = 1 To N
            PCOM(J) = P(J)
            XICOM(J) = XI(J)
        Next J
        AX = 0#
        XX = 1#
        Call MNBRAK(AX, XX, BX, FA, FX, FB)
        FRET = BRENT(AX, XX, BX, TOL, XMIN)
        For J = 1 To N
            XI(J) = XMIN * XI(J)
            P(J) = P(J) + XI(J)
        Next J
    End Sub
    Function F1DIM(X)
        Dim XT(50)
        For J = 1 To NCOM
            XT(J) = PCOM(J) + X * XICOM(J)
        Next J
        F1DIM = FUNC2(XT(), NCOM)
        Erase XT
    End Function
Sub MNBRAK(AX, BX, CX, FA, FB, FC)
    GOLD = 1.618034
    GLIMIT = 100#
    TINY = 1E-20
    FA = FUNC(AX)
    FB = FUNC(BX)
    If FB > FA Then
        DUM = AX
        AX = BX
        BX = DUM
        DUM = FB
        FB = FA
        FA = DUM
    End If
    CX = BX + GOLD * (BX - AX)
    FC = FUNC(CX)
    Do
        If FB < FC Then Exit Do
        DONE% = -1
        R = (BX - AX) * (FB - FC)
        Q = (BX - CX) * (FB - FA)
        DUM = Q - R
        If Abs(DUM) < TINY Then DUM = TINY
        U = BX - ((BX - CX) * Q - (BX - AX) * R) / (2 * DUM)
        ULIM = BX + GLIMIT * (CX - BX)
        If (BX - U) * (U - CX) > 0# Then
            FU = FUNC(U)
            If FU < FC Then
                AX = BX
                FA = FB
                BX = U
                FB = FU
                Exit Sub
            ElseIf FU > FB Then
                CX = U
                FC = FU
                Exit Sub
            End If
            U = CX + GOLD * (CX - BX)
            FU = FUNC(U)
        ElseIf (CX - U) * (U - ULIM) > 0# Then
            FU = FUNC(U)
            If FU < FC Then
                BX = CX
                CX = U
                U = CX + GOLD * (CX - BX)
                FB = FC
                FC = FU
                FU = FUNC(U)
            End If
        ElseIf (U - ULIM) * (ULIM - CX) >= 0# Then
            U = ULIM
            FU = FUNC(U)
        Else
            U = CX + GOLD * (CX - BX)
            FU = FUNC(U)
        End If
        If DONE% Then
            AX = BX
            BX = CX
            CX = U
            FA = FB
            FB = FC
            FC = FU
        Else
            DONE% = 0
        End If
    Loop While Not DONE%
End Sub
Function BRENT(AX, BX, CX, TOL, XMIN)
    ITMAX = 100
    CGOLD = 0.381966
    ZEPS = 0.0000000001
    A = AX
    If CX < AX Then A = CX
    B = AX
    If CX > AX Then B = CX
    V = BX
    W = V
    X = V
    E = 0#
    FX = FUNC(X)
    FV1 = FX
    FW = FX
    For ITER = 1 To ITMAX
        XM = 0.5 * (A + B)
        TOL1 = TOL * Abs(X) + ZEPS
        TOL2 = 2# * TOL1
        If Abs(X - XM) <= TOL2 - 0.5 * (B - A) Then Exit For
        DONE% = -1
        If Abs(E) > TOL1 Then
            R = (X - W) * (FX - FV1)
            Q = (X - V) * (FX - FW)
            P = (X - V) * Q - (X - W) * R
            Q = 2# * (Q - R)
            If Q > 0# Then P = -P
            Q = Abs(Q)
            ETEMP = E
            E = D
            DUM = Abs(0.5 * Q * ETEMP)
            If Abs(P) < DUM And P > Q * (A - X) And P < Q * (B - X) Then
                D = P / Q
                U = X + D
                If U - A < TOL2 Or B - U < TOL2 Then
                    D = Abs(TOL1) * Sgn(XM - X)
                End If
                DONE% = 0
            End If
        End If
        If DONE% Then
            If X >= XM Then
                E = A - X
            Else
                E = B - X
            End If
            D = CGOLD * E
        End If
        If Abs(D) >= TOL1 Then
            U = X + D
        Else
            U = X + Abs(TOL1) * Sgn(D)
        End If
        FU = FUNC(U)
        If FU <= FX Then
            If U >= X Then
                A = X
            Else
                B = X
            End If
            V = W
            FV1 = FW
            W = X
            FW = FX
            X = U
            FX = FU
        Else
            If U < X Then
                A = U
            Else
                B = U
            End If
            If FU <= FW Or W = X Then
                V = W
                FV1 = FW
                W = U
                FW = FU
            ElseIf FU <= FV1 Or V = X Or V = W Then
                V = U
                FV1 = FU
            End If
        End If
    Next ITER
    If ITER > ITMAX Then Form1.Print " Brent exceed maximum iterations."
    XMIN = X
    BRENT = FX
End Function
Function FUNC2(X(), N)
    FUNC2 = 1# - BESSJ0(X(1) - 0.5) * BESSJ0(X(2) - 0.5) * BESSJ0(X(3) - 0.5)
End Function
Sub DFUNC(X(), DF())
    NMAX = 50
    DF(1) = BESSJ1(X(1) - 0.5) * BESSJ0(X(2) - 0.5) * BESSJ0(X(3) - 0.5)
    DF(2) = BESSJ0(X(1) - 0.5) * BESSJ1(X(2) - 0.5) * BESSJ0(X(3) - 0.5)
    DF(3) = BESSJ0(X(1) - 0.5) * BESSJ0(X(2) - 0.5) * BESSJ1(X(3) - 0.5)
End Sub
Function FUNC(X)
    FUNC = F1DIM(X)
End Function
Function BESSJ0(X)
    P1# = 1#
    P2# = -0.001098628627
    P3# = 0.00002734510407
    P4# = -0.000002073370639
    P5# = 2.093887211E-07
    Q1# = -0.01562499995
    Q2# = 0.0001430488765
    Q3# = -0.000006911147651
    Q4# = 7.621095161E-07
    Q5# = -9.34945152E-08
    R1# = 57568490574#
    R2# = -13362590354#
    R3# = 651619640.7
    R4# = -11214424.18
    R5# = 77392.33017
    R6# = -184.9052456
    S1# = 57568490411#
    S2# = 1029532985#
    S3# = 9494680.718
    S4# = 59272.64853
    S5# = 267.8532712
    S6# = 1#
    If Abs(X) < 8# Then
       Y# = X * X
       BBB# = Y# * (R4# + Y# * (R5# + Y# * R6#))
       AAA# = R1# + Y# * (R2# + Y# * (R3# + BBB#))
       CCC# = Y# * (S3# + Y# * (S4# + Y# * (S5# + Y# * S6#)))
       BESSJ0 = AAA / (S1# + Y# * (S2# + CCC#))
    Else
       AX# = Abs(X)
       Z# = 8# / AX#
       Y# = Z# * Z#
       XX = AX# - 0.785398164
       CCC# = Y# * (P3# + Y# * (P4# + Y# * P5#))
       AAA# = P1# + Y# * (P2# + CCC#)
       DDD# = Y# * (Q3# + Y# * (Q4# + Y# * Q5#))
       EEE# = Z# * Sin(XX) * (Q1# + Y# * (Q2# + DDD#))
       BESSJ0 = Sqr(0.636619772 / AX#) * (Cos(XX) * AAA# - EEE#)
    End If
End Function
   Function BESSJ1(X)
      R1 = 72362614232#
      R2 = -7895059235#
      R3 = 242396853.1
      R4 = -2972611.439
      R5 = 15704.4826
      R6 = -30.16036606
      S1 = 144725228442#
      S2 = 2300535178#
      S3 = 18583304.74
      S4 = 99447.43394
      S5 = 376.9991397
      S6 = 1#
      P1 = 1#
      P2 = 0.00183105
      P3 = -0.00003516396496
      P4 = 0.000002457520174
      P5 = -0.000000240337019
      Q1 = 0.04687499995
      Q2 = -0.0002002690873
      Q3 = 0.000008449199096
      Q4 = -0.00000088228987
      Q5 = 0.000000105787412
      If Abs(X) < 8# Then
          Y# = X ^ 2
          AAA# = R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))))
          BBB# = S1 + Y * (S2 + Y * (S3 + Y * (S4 + Y * (S5 + Y * S6))))
          BESSJ1 = X * AAA# / BBB#
      Else
          AX = Abs(X)
          Z = 8# / AX
          Y# = Z ^ 2
          XX = AX - 2.356194491
          AAA# = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
          BBB# = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
          CCC# = Sqr(0.636619772 / AX)
          BESSJ1 = CCC# * (Cos(XX) * AAA# - Z * Sin(XX) * BBB# * Sgn(X))
      End If
   End Function