d3r4.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3195
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   3195
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      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 D3R4
    'Driver for routine QROMB
    PIO2 = 1.5707963
    A = 0#
    B = PIO2
    Print
    Print Tab(5); "Integral of FUNC with QROMB"
    Print
    Print Tab(5); "Actual value of integral is:";
    Print Tab(34); Format$((FINT(B) - FINT(A)), "0.####00")
    Print
    Call QROMB(A, B, S)
    Print Tab(5); "Result from routine QROMB is:";
    Print Tab(35); Format$(S, "0.####00")
End Sub
Function FUNC(X)
    FUNC = (X ^ 2) * (X ^ 2 - 2#) * Sin(X)
End Function
Function FINT(X)
    'Integral of FUNC
    AAA = 4# * X * (X ^ 2 - 7#) * Sin(X)
    FINT = AAA - (X ^ 4 - 14# * X ^ 2 + 28#) * Cos(X)
End Function
Sub QROMB(A, B, SS)
    EPS = 0.000001
    JMAX = 20
    JMAXP = JMAX + 1
    K = 5
    KM = K - 1
    Dim S(21), H(21)
    H(1) = 1#
    For J = 1 To JMAX
        Call TRAPZD(A, B, S(J), J)
        If J >= K Then
            Call POLINT(H(), S(), K, 0#, SS, DSS)
            If Abs(DSS) < EPS * Abs(SS) Then Exit Sub
        End If
        S(J + 1) = S(J)
        H(J + 1) = 0.25 * H(J)
    Next J
    Print "Too many steps."
End Sub
    Sub TRAPZD(A, B, S, N)
       If N = 1 Then
           S = 0.5 * (B - A) * (FUNC(A) + FUNC(B))
           IT = 1
       Else
           IT = 2 ^ (N - 2)
           TNM = IT
           DEL = (B - A) / TNM
           X = A + 0.5 * DEL
           Sum = 0#
           For J = 1 To IT
               Sum = Sum + FUNC(X)
               X = X + DEL
           Next J
           S = 0.5 * (S + (B - A) * Sum / TNM)
       End If
    End Sub
   Sub POLINT(XA(), YA(), N, X, Y, DY)
      Dim C(10), D(10)
      NS = 1
      DIF = Abs(X - XA(1))
      For I = 1 To N
          DIFT = Abs(X - XA(I))
          If DIFT < DIF Then
              NS = I
              DIF = DIFT
          End If
          C(I) = YA(I)
          D(I) = YA(I)
      Next I
      Y = YA(NS)
      NS = NS - 1
      For M = 1 To N - 1
          For I = 1 To N - M
              HO = XA(I) - X
              HP = XA(I + M) - X
              W = C(I + 1) - D(I)
              DEN = HO - HP
              If DEN = 0# Then
                 Print "PAUSE"
                 Exit Sub
              End If
              DEN = W / DEN
              D(I) = HP * DEN
              C(I) = HO * DEN
          Next I
          If 2 * NS < N - M Then
              DY = C(NS + 1)
          Else
              DY = D(NS)
              NS = NS - 1
          End If
          Y = Y + DY
      Next M
   End Sub