d3r9.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3555
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   3555
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3000
      TabIndex        =   0
      Top             =   3000
      Width           =   1335
   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 D3R9
    'Driver for routine GAULEG
    NPOINT = 10
    X1 = 0#
    X2 = 1#
    X3 = 10#
    Dim X(10), W(10)
    Call GAULEG(X1, X2, X(), W(), NPOINT)
    Print
    Print Tab(5); "#           X(I)        W(I)"
    For I = 1 To NPOINT
        Print Tab(5); Format$(I, "##");
        Print Tab(14); Format$(X(I), "0.######");
        Print Tab(27); Format$(W(I), "0.######")
    Next I
    'Demonstrate the use of GAUEG for an integral
    Call GAULEG(X1, X3, X(), W(), NPOINT)
    XX = 0#
    For I = 1 To NPOINT
        XX = XX + W(I) * FUNC(X(I))
    Next I
    Print
    Print Tab(5); "Integral from GAULEG: "; Format$(XX, "0.######")
    Print Tab(5); "Actual value: ";
    Print Tab(20); Format$((1# - (1# + X3) * Exp(-X3)), "0.######")
End Sub
Function FUNC(X)
    FUNC = X * Exp(-X)
End Function
Sub GAULEG(X1, X2, X(), W(), N)
    EPS = 0.00000000000003
    M = (N + 1) / 2
    XM = 0.5 * (X2 + X1)
    XL = 0.5 * (X2 - X1)
    For I = 1 To M
        Z = Cos(3.141592654 * (I - 0.25) / (N + 0.5))
1       P1 = 1#
        P2 = 0#
        For J = 1 To N
            P3 = P2
            P2 = P1
            P1 = ((2# * J - 1#) * Z * P2 - (J - 1#) * P3) / J
        Next J
        PP = N * (Z * P1 - P2) / (Z * Z - 1#)
        Z1 = Z
        Z = Z1 - P1 / PP
        If Abs(Z - Z1) > EPS Then GoTo 1
        X(I) = XM - XL * Z
        X(N + 1 - I) = XM + XL * Z
        W(I) = 2# * XL / ((1# - Z * Z) * PP * PP)
        W(N + 1 - I) = W(I)
    Next I
End Sub