d1r7.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   4875
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   4875
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   2880
      TabIndex        =   0
      Top             =   4320
      Width           =   1575
   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 D1R7
    'Driver program for routine TOEPLZ
    N = 5
    N2 = 2 * N
    Dim X(5), Y(5), R(10)
    '输入已知的方程组的系数矩阵的R
    For I = 1 To 2 * N - 1
    R(I) = 0.1 * I
    Next I
    '输入已知的方程组的右端向量
    For I = 1 To N
        Y(I) = 0.1 * I
    Next I
    Print Tab(5); "已知的方程组的右端向量"
    Print Tab(12); Format$(Y(1), "##.##")
    Print Tab(12); Format$(Y(2), "##.##")
    Print Tab(12); Format$(Y(3), "##.##")
    Print Tab(12); Format$(Y(4), "##.##")
    Print Tab(12); Format$(Y(5), "##.##")
    Call TOEPLZ(R(), X(), Y(), N)
    '输出方程组的解X
    Print
    Print Tab(5); "计算出的方程组的的解"
    Print Tab(9); Format$(X(1), "##.######")
    Print Tab(9); Format$(X(2), "##.######")
    Print Tab(9); Format$(X(3), "##.######")
    Print Tab(9); Format$(X(4), "##.######")
    Print Tab(9); Format$(X(5), "##.######")
    '将计算出的解X乘以系数矩阵，以验证计算结果正确
    Print
    Print Tab(4); "将计算出的解乘以系数矩阵，以验证计算结果正确"
    Print
    Print Tab(5); "解乘以系数矩阵"; Tab(24); "方程组的右端向量"
    For I = 1 To N
        SUM1 = 0
        For J = 1 To N
           SUM1 = SUM1 + R(N + I - J) * X(J)
        Next J
        Print Tab(10); Format$(SUM1, "#.####"); Tab(30); Format$(Y(I), "#.####")
    Next I
End Sub
Sub TOEPLZ(R(), X(), Y(), N)
    NMAX = 100
    Dim G(100), H(100)
    If R(N) = 0# Then GoTo 99
    X(1) = Y(1) / R(N)
    If N = 1 Then Exit Sub
    G(1) = R(N - 1) / R(N)
    H(1) = R(N + 1) / R(N)
    For M = 1 To N
        M1 = M + 1
        SXN = -Y(M1)
        SD = -R(N)
        For J = 1 To M
            SXN = SXN + R(N + M1 - J) * X(J)
            SD = SD + R(N + M1 - J) * G(M - J + 1)
        Next J
        If SD = 0# Then GoTo 99
        X(M1) = SXN / SD
        For J = 1 To M
            X(J) = X(J) - X(M1) * G(M - J + 1)
        Next J
        If M1 = N Then Exit Sub
        SGN1 = -R(N - M1)
        SHN = -R(N + M1)
        SGD = -R(N)
        For J = 1 To M
            SGN1 = SGN1 + R(N + J - M1) * G(J)
            SHN = SHN + R(N + M1 - J) * H(J)
            SGD = SGD + R(N + J - M1) * H(M - J + 1)
        Next J
        If SD = 0# Or SGD = 0# Then GoTo 99
        G(M1) = SGN1 / SGD
        H(M1) = SHN / SD
        K = M
        M2 = (M + 1) / 2
        PP = G(M1)
        QQ = H(M1)
        For J = 1 To M2
            PT1 = G(J)
            PT2 = G(K)
            QT1 = H(J)
            QT2 = H(K)
            G(J) = PT1 - PP * QT2
            G(K) = PT2 - PP * QT1
            H(J) = QT1 - QQ * PT2
            H(K) = QT2 - QQ * PT1
            K = K - 1
        Next J
    Next M
    Print "Never get here"
99  Print "Levinson method fails: singular principal minor"
End Sub