d9r5.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3195
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   5550
   LinkTopic       =   "Form1"
   ScaleHeight     =   3195
   ScaleWidth      =   5550
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3720
      TabIndex        =   0
      Top             =   2640
      Width           =   1095
   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 D9R5
    'Driver for routine AMOEBA
    NP = 3
    MP = 4
    FTOL = 0.000001
    Dim P(4, 3), X(3), Y(4)
    P(1, 1) = 0#: P(1, 2) = 0#: P(1, 3) = 0#
    P(2, 1) = 1#: P(2, 2) = 0#: P(2, 3) = 0#
    P(3, 1) = 0#: P(3, 2) = 1#: P(3, 3) = 0#
    P(4, 1) = 0#: P(4, 2) = 0#: P(4, 3) = 1#
    NDIM = NP
    For I = 1 To MP
        For J = 1 To NP
            X(J) = P(I, J)
        Next J
        Y(I) = FAMOEB(X)
    Next I
    Call AMOEBA(P(), Y(), MP, NP, NDIM, FTOL, ITER)
    Print
    Print Tab(5); "Iterations: ", Format$(ITER, "##")
    Print Tab(5); "Vertices of final 3-D simplex and"
    Print Tab(5); "function values at the vertices:"
    Print Tab(5); "I        X(I)       Y(I)       Z(I)      FUNCTION"
    For I = 1 To MP
        Print Tab(5); Format$(I, "#");
        For J = 1 To NP
            Print Tab(J * 11); Format$(P(I, J), ".#####0");
        Next J
        Print Tab(46); Format$(Y(I), ".#####0")
    Next I
    Print
    Print Tab(5); "True minimum is at ( 0.5, 0.6, 0.7)"
End Sub
Function FAMOEB(X())
    AAA = (X(1) - 0.5) ^ 2 + (X(2) - 0.6) ^ 2 + (X(3) - 0.7) ^ 2
    FAMOEB = 0.6 - BESSJ0(AAA)
End Function
Sub AMOEBA(P(), Y(), MP, NP, NDIM, FTOL, ITER)
    NMAX = 20
    ALPHA = 1#
    BETA = 0.5
    GAMMA = 2#
    ITMAX = 500
    Dim PR(20), PRR(20), PBAR(20)
    MPTS = NDIM + 1
    ITER = 0
    Do
      ILO = 1
      If Y(1) > Y(2) Then
          IHI = 1
          INHI = 2
      Else
          IHI = 2
          INHI = 1
      End If
      For I = 1 To MPTS
          If Y(I) < Y(ILO) Then ILO = I
          If Y(I) > Y(IHI) Then
              INHI = IHI
              IHI = I
          ElseIf Y(I) > Y(INHI) Then
              If I <> IHI Then INHI = I
          End If
      Next I
      RTOL = 2# * Abs(Y(IHI) - Y(ILO)) / (Abs(Y(IHI)) + Abs(Y(ILO)))
      If RTOL < FTOL Then Erase PBAR, PRR, PR: Exit Sub
      If ITER = ITMAX Then
          Print " Amoeba exceeding maximum iterations."
          Exit Sub
      End If
      ITER = ITER + 1
      For J = 1 To NDIM
          PBAR(J) = 0#
      Next J
      For I = 1 To MPTS
          If I <> IHI Then
              For J = 1 To NDIM
                  PBAR(J) = PBAR(J) + P(I, J)
              Next J
          End If
      Next I
      For J = 1 To NDIM
          PBAR(J) = PBAR(J) / NDIM
          PR(J) = (1# + ALPHA) * PBAR(J) - ALPHA * P(IHI, J)
      Next J
      YPR = FAMOEB(PR)
      If YPR <= Y(ILO) Then
          For J = 1 To NDIM
              PRR(J) = GAMMA * PR(J) + (1# - GAMMA) * PBAR(J)
          Next J
          YPRR = FAMOEB(PRR)
          If YPRR < Y(ILO) Then
              For J = 1 To NDIM
                  P(IHI, J) = PRR(J)
              Next J
              Y(IHI) = YPRR
          Else
              For J = 1 To NDIM
                  P(IHI, J) = PR(J)
              Next J
              Y(IHI) = YPR
          End If
      ElseIf YPR >= Y(INHI) Then
          If YPR < Y(IHI) Then
              For J = 1 To NDIM
                  P(IHI, J) = PR(J)
              Next J
              Y(IHI) = YPR
          End If
          For J = 1 To NDIM
              PRR(J) = BETA * P(IHI, J) + (1# - BETA) * PBAR(J)
          Next J
          YPRR = FAMOEB(PRR)
          If YPRR < Y(IHI) Then
              For J = 1 To NDIM
                  P(IHI, J) = PRR(J)
              Next J
              Y(IHI) = YPRR
          Else
              For I = 1 To MPTS
                  If I <> ILO Then
                      For J = 1 To NDIM
                          PR(J) = 0.5 * (P(I, J) + P(ILO, J))
                          P(I, J) = PR(J)
                      Next J
                      Y(I) = FAMOEB(PR)
                  End If
              Next I
          End If
      Else
          For J = 1 To NDIM
              P(IHI, J) = PR(J)
          Next J
          Y(IHI) = YPR
      End If
    Loop
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