d9r1.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6615
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   5415
   LinkTopic       =   "Form1"
   ScaleHeight     =   6615
   ScaleWidth      =   5415
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3960
      TabIndex        =   0
      Top             =   5880
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
    'PROGRAM D9R1
    'Driver for routine MNBRAK
    For I = 1 To 10
        AX = I * 0.5
        BX = (I + 1#) * 0.5
        Call MNBRAK(AX, BX, CX, FA, FB, FC)
        Print Tab(5); "              A             B             C"
        Print Tab(5); "X";
        Print Tab(15); Format$(AX, "0.####00");
        Print Tab(29); Format$(BX, "0.####00");
        Print Tab(43); Format$(CX, "0.####00")
        Print Tab(5); "F";
        Print Tab(15); Format$(FA, "0.####00");
        Print Tab(29); Format$(FB, "0.####00");
        Print Tab(43); Format$(FC, "0.####00")
    Next I
End Sub
Function FUNC(X)
    FUNC = BESSJ0(X)
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 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