d3r6.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   4770
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   5295
   LinkTopic       =   "Form1"
   ScaleHeight     =   4770
   ScaleWidth      =   5295
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3600
      TabIndex        =   0
      Top             =   4080
      Width           =   1335
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Public CHOOSE As String
Private Sub Command1_Click()
    'PROGRAM D3R6
    'Driver for routine QROMO
    X1 = 0#
    X2 = 1.5707963
    X3 = 3.1415926
    AINF = 1E+20
    Print
    Print Tab(5); "Improper integrals:"
    Print
    CHOOSE = "FUNCL"
    Call QROMO(X1, X2, RESULT, "MIDSQL")
    Print Tab(5); "Function: SQR(x)/SIN(x)      Interval: (0,PI/2)"
    Print Tab(5); "Using: MIDSQL                Result: ";
    Print Format$(RESULT, "0.####")
    Print
    CHOOSE = "FUNCU"
    Call QROMO(X2, X3, RESULT, "MIDSQU")
    Print Tab(5); "Function: SQR(PI-x)/SIN(x)   Interval: (PI/2,PI)"
    Print Tab(5); "Using: MIDSQU                Result: ";
    Print Format$(RESULT, "0.####")
    Print
    CHOOSE = "FUNCINF"
    Call QROMO(X2, AINF, RESULT, "MIDINF")
    Print Tab(5); "Function: SIN(x)/x^2        Interval: (PI/2,infty)"
    Print Tab(5); "Using: MIDINF                Result: ";
    Print Format$(RESULT, "0.####")
    Print
    CHOOSE = "FUNCINF"
    Call QROMO(-AINF, -X2, RESULT, "MIDINF")
    Print Tab(5); "Function: SIN(x)/x^2        Interval: (-infty,-PI/2)"
    Print Tab(5); "Using: MIDINF                Result: ";
    Print Format$(RESULT, "0.####")
    Print
    CHOOSE = "FUNCEND"
    Call QROMO(X1, X2, RES1, "MIDSQL")
    Call QROMO(X2, AINF, RES2, "MIDINF")
    Print Tab(5); "Function: EXP(-x)/SQR(x)     Interval: (0,infty)"
    Print Tab(5); "Using: MIDSQL,MIDINF         Result: ";
    Print Format$(RES1 + RES2, "0.####")
End Sub
Function FUNC(X)
    'Dim CHOOSE As String
    If CHOOSE = "FUNCL" Then FUNC = FUNCL(X)
    If CHOOSE = "FUNCU" Then FUNC = FUNCU(X)
    If CHOOSE = "FUNCINF" Then FUNC = FUNCINF(X)
    If CHOOSE = "FUNCEND" Then FUNC = FUNCEND(X)
End Function
Function FUNCL(X)
    FUNCL = Sqr(X) / Sin(X)
End Function
Function FUNCU(X)
    PI = 3.1415926
    FUNCU = Sqr(PI - X) / Sin(X)
End Function
Function FUNCINF(X)
    FUNCINF = Sin(X) / (X ^ 2)
End Function
Function FUNCEND(X)
    FUNCEND = Exp(-X) / Sqr(X)
End Function
Sub QROMO(A, B, SS, PICK$)
    EPS = 0.00003
    JMAX = 14
    JMAXP = JMAX + 1
    K = 5
    KM = K - 1
    Dim S(15), H(15)
    H(1) = 1#
    For J = 1 To JMAX
        If PICK$ = "MIDPNT" Then Call MIDPNT(A, B, S(J), J)
        If PICK$ = "MIDINF" Then Call MIDINF(A, B, S(J), J)
        If PICK$ = "MIDSQL" Then Call MIDSQL(A, B, S(J), J)
        If PICK$ = "MIDSQU" Then Call MIDSQU(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) = H(J) / 9#
    Next J
    Print "Too many steps."
End Sub
Sub MIDINF(AA, BB, S, N)
    B = 1# / AA
    A = 1# / BB
    If N = 1 Then
        X = 0.5 * (A + B)
        S = (B - A) * FUNC(0.5 * (A + B))
        IT = 1
    Else
        IT = 3 ^ (N - 2)
        TNM = IT
        DEL = (B - A) / (3# * TNM)
        DDEL = DEL + DEL
        X = A + 0.5 * DEL
        Sum = 0#
        For J = 1 To IT
            Sum = Sum + INF(X)
            X = X + DDEL
            Sum = Sum + INF(X)
            X = X + DEL
        Next J
        S = (S + (B - A) * Sum / TNM) / 3#
    End If
End Sub
Function INF(X)
    INF = FUNC(1 / X) / X ^ 2
End Function
Sub MIDSQL(AA, BB, S, N)
    B = Sqr(BB - AA)
    A = 0#
    If N = 1 Then
        S = (B - A) * SQL(0.5 * (A + B), AA)
        IT = 1
    Else
        IT = 3 ^ (N - 2)
        TNM = IT
        DEL = (B - A) / (3# * TNM)
        DDEL = DEL + DEL
        X = A + 0.5 * DEL
        Sum = 0#
        For J = 1 To IT
            Sum = Sum + SQL(X, AA)
            X = X + DDEL
            Sum = Sum + SQL(X, AA)
            X = X + DEL
        Next J
        S = (S + (B - A) * Sum / TNM) / 3#
    End If
End Sub
Function SQL(X, AA)
    SQL = 2 * X * FUNC(AA + X ^ 2)
End Function
Sub MIDSQU(AA, BB, S, N)
    B = Sqr(BB - AA)
    A = 0#
    If N = 1 Then
        S = (B - A) * SQU(0.5 * (A + B), BB)
        IT = 1
    Else
        IT = 3 ^ (N - 2)
        TNM = IT
        DEL = (B - A) / (3# * TNM)
        DDEL = DEL + DEL
        X = A + 0.5 * DEL
        Sum = 0#
        For J = 1 To IT
            Sum = Sum + SQU(X, BB)
            X = X + DDEL
            Sum = Sum + SQU(X, BB)
            X = X + DEL
        Next J
        S = (S + (B - A) * Sum / TNM) / 3#
    End If
End Sub
Function SQU(X, BB)
    SQU = 2 * X * FUNC(BB - X ^ 2)
End Function
Sub MIDPNT(A, B, S, N)
    If N = 1 Then
        S = (B - A) * FUNC(0.5 * (A + B))
        IT = 1
    Else
        IT = 3 ^ (N - 2)
        TNM = IT
        DEL = (B - A) / (3# * TNM)
        DDEL = DEL + DEL
        X = A + 0.5 * DEL
        Sum = 0#
        For J = 1 To IT
            Sum = Sum + FUNC(X)
            X = X + DDEL
            Sum = Sum + FUNC(X)
            X = X + DEL
        Next J
        S = (S + (B - A) * Sum / TNM) / 3#
    End If
End Sub
Sub POLINT(XA(), YA(), N, X, Y, DY)
    Dim C(15), D(15)
    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