d1r5.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   4620
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   4620
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3120
      TabIndex        =   0
      Top             =   3960
      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 D1R5
    'Driver program for routine MPROVE
    N = 5
    Dim A(5, 5), A1(5, 5), B(5), INDX(5), B1(5)
    '输入已知的方程组的系数矩阵
    A(1, 1) = 1: A(1, 2) = 2: A(1, 3) = 3: A(1, 4) = 4: A(1, 5) = 5
    A(2, 1) = 2: A(2, 2) = 3: A(2, 3) = 4: A(2, 4) = 5: A(2, 5) = 1
    A(3, 1) = 1: A(3, 2) = 1: A(3, 3) = 1: A(3, 4) = 1: A(3, 5) = 1
    A(4, 1) = 4: A(4, 2) = 5: A(4, 3) = 1: A(4, 4) = 2: A(4, 5) = 3
    A(5, 1) = 5: A(5, 2) = 1: A(5, 3) = 2: A(5, 4) = 3: A(5, 5) = 4
    '输入已知的方程组的右端向量B
    B(1) = 1
    B(2) = 1
    B(3) = 1
    B(4) = 1
    B(5) = 1
    For I = 1 To N
        For J = 1 To N
            A1(I, J) = A(I, J)
        Next J
    Next I
    For I = 1 To N
        B1(I) = B(I)
    Next I
    Call LUDCMP(A(), N, INDX(), D)
    Call LUBKSB(A(), N, INDX(), B())
    '输出方程组的解B
    Print
    Print Tab(5); "输出方程组的解"
    Print Tab(10); Format$(B(1), "##.####")
    Print Tab(10); Format$(B(2), "##.####")
    Print Tab(10); Format$(B(3), "##.####")
    Print Tab(10); Format$(B(4), "##.####")
    Print Tab(10); Format$(B(5), "##.####")
    For I = 1 To N
        B(I) = B(I) * (1# + 0.2 * Rnd(I))
    Next I
    '输出干扰后的解B
    Print
    Print Tab(5); "输出干扰后的解"
    Print Tab(8); Format$(B(1), "##.####")
    Print Tab(8); Format$(B(2), "##.####")
    Print Tab(8); Format$(B(3), "##.####")
    Print Tab(8); Format$(B(4), "##.####")
    Print Tab(8); Format$(B(5), "##.####")
    Call MPROVE(A1(), A(), N, INDX(), B1(), B())
    '输出改善后的解B
    Print
    Print Tab(5); "输出改善后的解"
    Print Tab(10); Format$(B(1), "##.####")
    Print Tab(10); Format$(B(2), "##.####")
    Print Tab(10); Format$(B(3), "##.####")
    Print Tab(10); Format$(B(4), "##.####")
    Print Tab(10); Format$(B(5), "##.####")
End Sub
Sub MPROVE(A(), ALUD(), N, INDX(), B(), X())
    NMAX = 100
    Dim R(100)
    For I = 1 To N
        SDP = -B(I)
        For J = 1 To N
            SDP = SDP + A(I, J) * X(J)
        Next J
        R(I) = SDP
    Next I
    Call LUBKSB(ALUD(), N, INDX(), R())
    For I = 1 To N
        X(I) = X(I) - R(I)
    Next I
End Sub
   Sub LUBKSB(A(), N, INDX, B())
      II = 0
      For I = 1 To N
          LL = INDX(I)
          Sum = B(LL)
          B(LL) = B(I)
          If II <> 0 Then
              For J = II To I - 1
                  Sum = Sum - A(I, J) * B(J)
              Next J
          ElseIf Sum <> 0# Then
              II = I
          End If
          B(I) = Sum
      Next I
      For I = N To 1 Step -1
          Sum = B(I)
          If I < N Then
              For J = I + 1 To N
                  Sum = Sum - A(I, J) * B(J)
              Next J
          End If
          B(I) = Sum / A(I, I)
      Next I
   End Sub
   Sub LUDCMP(A(), N, INDX(), D)
      NMAX = 100
      TINY = 1E-20
      Dim VV(100)
      D = 1#
      For I = 1 To N
          AAMAX = 0#
          For J = 1 To N
              If Abs(A(I, J)) > AAMAX Then AAMAX = Abs(A(I, J))
          Next J
          If AAMAX = 0# Then Print "Singular matrix."
          VV(I) = 1# / AAMAX
      Next I
      For J = 1 To N
          If J > 1 Then
              For I = 1 To J - 1
                  Sum = A(I, J)
                  If I > 1 Then
                      For K = 1 To I - 1
                          Sum = Sum - A(I, K) * A(K, J)
                      Next K
                      A(I, J) = Sum
                  End If
              Next I
          End If
          AAMAX = 0#
          For I = J To N
              Sum = A(I, J)
              If J > 1 Then
                  For K = 1 To J - 1
                      Sum = Sum - A(I, K) * A(K, J)
                  Next K
                  A(I, J) = Sum
              End If
              DUM = VV(I) * Abs(Sum)
              If DUM >= AAMAX Then
                  IMAX = I
                  AAMAX = DUM
              End If
          Next I
          If J <> IMAX Then
              For K = 1 To N
                  DUM = A(IMAX, K)
                  A(IMAX, K) = A(J, K)
                  A(J, K) = DUM
              Next K
              D = -D
              VV(IMAX) = VV(J)
          End If
          INDX(J) = IMAX
          If J <> N Then
              If A(J, J) = 0# Then A(J, J) = TINY
              DUM = 1# / A(J, J)
              For I = J + 1 To N
                  A(I, J) = A(I, J) * DUM
              Next I
          End If
      Next J
      If A(N, N) = 0# Then A(N, N) = TINY
   End Sub