d11r11.frm

方程组的矩阵解法

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3765
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   3765
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3120
      TabIndex        =   0
      Top             =   3000
      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 D11R11
    'Driver for routine KSTWO
    N1 = 2000
    N2 = 1000
    EPS = 0.1
    Dim DATA1(2000), DATA2(1000)
    IDUM = -1357
    For J = 1 To N1
        DATA1(J) = GASDEV(IDUM)
    Next J
    Print
    Print Tab(5); "Variance Ratio    K-S Statistic   Probability"
    IDUM = -2468
    For I = 1 To 11
        Var = 1# + (I - 1) * EPS
        FACTR = Sqr(Var)
        For J = 1 To N2
            DATA2(J) = FACTR * GASDEV(IDUM)
        Next J
        Call KSTWO(DATA1(), N1, DATA2(), N2, D, PROB)
        Print Tab(8); Format$(Var, "#.#####0");
        Print Tab(23); Format$(D, "#.#####0");
        Print Tab(39); Format$(PROB, ".####E+00")
    Next I
End Sub
Sub KSTWO(DATA1(), N1, DATA2(), N2, D, PROB)
    Call SORT(N1, DATA1)
    Call SORT(N2, DATA2)
    EN1 = N1
    EN2 = N2
    J1 = 1
    J2 = 1
    FO1 = 0#
    FO2 = 0#
    D = 0#
1   If J1 <= N1 And J2 <= N2 Then
        If DATA1(J1) < DATA2(J2) Then
            FN1 = J1 / EN1
            If Abs(FN1 - FO2) > Abs(FO1 - FO2) Then
                AM = Abs(FN1 - FO2)
            Else
                AM = Abs(FO1 - FO2)
            End If
            DT = AM
            If DT > D Then D = DT
            FO1 = FN1
            J1 = J1 + 1
        Else
            FN2 = J2 / EN2
            If Abs(FN2 - FO1) > Abs(FO2 - FO1) Then
                AM = Abs(FN2 - FO1)
            Else
                AM = Abs(FO2 - FO1)
            End If
            DT = AM
            If DT > D Then D = DT
            FO2 = FN2
            J2 = J2 + 1
        End If
        GoTo 1
    End If
    PROB = PROBKS(Sqr(EN1 * EN2 / (EN1 + EN2)) * D)
End Sub
   Function GASDEV(IDUM)
      Static ISET, GSET
      If ISET = 0 Then
        Do
          Randomize
          V1 = 2# * Rnd - 1#
          V2 = 2# * Rnd - 1#
          R = V1 ^ 2 + V2 ^ 2
        Loop While R >= 1# Or R = 0
        FAC = Sqr(-2# * Log(R) / R)
        GSET = V1 * FAC
        GASDEV = V2 * FAC
        ISET = 1
      Else
        GASDEV = GSET
        ISET = 0
      End If
   End Function
Sub SORT(N, RA())
      L = Int(N / 2) + 1
      IR = N
      Do
          If L > 1 Then
            L = L - 1
            RRA = RA(L)
          Else
            RRA = RA(IR)
            RA(IR) = RA(1)
            IR = IR - 1
            If IR = 1 Then
              RA(1) = RRA
              Exit Sub
            End If
          End If
          I = L
          J = L + L
          While J <= IR
            If J < IR Then
              If RA(J) < RA(J + 1) Then J = J + 1
            End If
            If RRA < RA(J) Then
              RA(I) = RA(J)
              I = J
              J = J + J
            Else
              J = IR + 1
            End If
          Wend
          RA(I) = RRA
      Loop
End Sub
Function PROBKS(ALAM)
      EPS1 = 0.001
      EPS2 = 0.00000001
      A2 = -2# * ALAM ^ 2
      FAC = 2#
      PROBKS = 0#
      TERMBF = 0#
      For J = 1 To 100
          TERM = FAC * Exp(A2 * J ^ 2)
          PROBKS = PROBKS + TERM
          AAA = Abs(TERM)
          If AAA < EPS1 * TERMBF Or AAA < EPS2 * PROBKS Then
              Exit Function
          End If
          FAC = -FAC
          TERMBF = Abs(TERM)
      Next J
      PROBKS = 1#
End Function