d9r9.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6210
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   6210
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      TabIndex        =   0
      Top             =   5520
      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 D9R9
      'Driver for routine DFPMIN
      NDIM = 3
      PIO2 = 1.5707963
      FTOL = 0.000001
      Dim P(3)
      Print Tab(5); "PROGRAM finds the minimum of a function"
      Print Tab(5); "with different trial starting vectors."
      Print Tab(5); "True minimum is (0.5, 0.5, 0.5)"
      For K = 0 To 4
          ANGL = PIO2 * K / 4#
          P(1) = 2# * Cos(ANGL)
          P(2) = 2# * Sin(ANGL)
          P(3) = 0#
          Print Tab(5)
          Print Tab(5); "Starting vector: ("; Format$(P(1), "#.#000");
          Print ","; Format$(P(2), "#.#000"); ","; Format$(P(3), "#.#000"); ")"
          Call DFPMIN(P(), NDIM, FTOL, ITER, FRET)
          Print Tab(5); "Iterations: "; Format$(ITER, "##")
          Print Tab(5); "Solution vector: ("; Format$(P(1), "#.#000");
          Print ","; Format$(P(2), "#.#000"); ","; Format$(P(3), "#.#000"); ")"
          Print Tab(5); "Func. value at solution", Format$(FRET, ".######E+00")
      Next K
End Sub
Sub DFUNC(X(), DF())
    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 FUNC2(X(), N)
    FUNC2 = 1# - BESSJ0(X(1) - 0.5) * BESSJ0(X(2) - 0.5) * BESSJ0(X(3) - 0.5)
End Function
Function FUNC(X)
    FUNC = F1DIM(X)
End Function
Sub DFPMIN(P(), N, FTOL, ITER, FRET)
    ITMAX = 200
    EPS = 0.0000000001
    Dim HESSIN(50, 50), XI(50), G(50), DG(50), HDG(50)
    FP = FUNC2(P(), N)
    Call DFUNC(P(), G())
    For I = 1 To N
        For J = 1 To N
            HESSIN(I, J) = 0#
        Next J
        HESSIN(I, I) = 1#
        XI(I) = -G(I)
    Next I
    For ITS = 1 To ITMAX
        ITER = ITS
        Call LINMIN(P(), XI(), N, FRET)
        If 2# * Abs(FRET - FP) <= FTOL * (Abs(FRET) + Abs(FP) + EPS) Then
            Erase HDG, DG, G, XI, HESSIN
            Exit Sub
        End If
        FP = FRET
        For I = 1 To N
            DG(I) = G(I)
        Next I
        FRET = FUNC2(P(), N)
        Call DFUNC(P(), G())
        For I = 1 To N
            DG(I) = G(I) - DG(I)
        Next I
        For I = 1 To N
            HDG(I) = 0#
            For J = 1 To N
                HDG(I) = HDG(I) + HESSIN(I, J) * DG(J)
            Next J
        Next I
        FAC = 0#
        FAE = 0#
        For I = 1 To N
            FAC = FAC + DG(I) * XI(I)
            FAE = FAE + DG(I) * HDG(I)
        Next I
        FAC = 1# / FAC
        FAD = 1# / FAE
        For I = 1 To N
            DG(I) = FAC * XI(I) - FAD * HDG(I)
        Next I
        For I = 1 To N
            For J = 1 To N
              AAA = FAC * XI(I) * XI(J) - FAD * HDG(I) * HDG(J)
              BBB = FAE * DG(I) * DG(J)
              HESSIN(I, J) = HESSIN(I, J) + AAA + BBB
            Next J
        Next I
        For I = 1 To N
            XI(I) = 0#
            For J = 1 To N
                XI(I) = XI(I) - HESSIN(I, J) * G(J)
            Next J
        Next I
    Next ITS
    If ITS > ITMAX Then Print " too many iterations in DFPMIN"
    Erase HDG, DG, G, XI, HESSIN
End Sub