d13r11.frm

VB常用数值算法集 内含有解线性代数方程组

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 D13R11
    '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
        V1 = 2# * RAN1(IDUM&) - 1#
        V2 = 2# * RAN1(IDUM&) - 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
Static Function RAN1(IDUM&)
    Dim R(97)
    M1& = 259200: IA1& = 7141: IC1& = 54773: RM1 = 0.0000038580247
    M2& = 134456: IA2& = 8121: IC2& = 28411: RM2 = 0.0000074373773
    M3& = 243000: IA3& = 4561: IC3& = 51349
    If IDUM& < 0 Or IFF = 0 Then
        IFF = 1
        IX1& = (IC1& - IDUM&) Mod M1&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX2& = IX1& Mod M2&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX3& = IX1& Mod M3&
        For J = 1 To 97
            IX1& = (IA1& * IX1& + IC1&) Mod M1&
            IX2& = (IA2& * IX2& + IC2&) Mod M2&
            R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
        Next J
        IDUM& = 1
    End If
    IX1& = (IA1& * IX1& + IC1&) Mod M1&
    IX2& = (IA2& * IX2& + IC2&) Mod M2&
    IX3& = (IA3& * IX3& + IC3&) Mod M3&
    J = 1 + Int((97 * IX3&) / M3&)
    If J > 97 Or J < 1 Then Print "Abnormal exit": Exit Function
    RAN1 = R(J)
    R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
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