d16r2.frm

VB常用数值算法集1 VB常用数值算法集1

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6270
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   8250
   LinkTopic       =   "Form1"
   ScaleHeight     =   6270
   ScaleWidth      =   8250
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   6000
      TabIndex        =   0
      Top             =   5280
      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 D16R2
    'Driver for routine ADI
    JMAX = 11
    PI = 3.1415926
    Dim A(11, 11), B(11, 11), C(11, 11), D(11, 11)
    Dim E(11, 11), F(11, 11), G(11, 11), U(11, 11)
    For I = 1 To JMAX
        For J = 1 To JMAX
            A(I, J) = -1#
            B(I, J) = 2#
            C(I, J) = -1#
            D(I, J) = -1#
            E(I, J) = 2#
            F(I, J) = -1#
            G(I, J) = 0#
            U(I, J) = 0#
        Next J
    Next I
    Mid1 = JMAX / 2 + 1
    G(Mid1, Mid1) = 2#
    ALPHA = 2# * (1# - Cos(PI / JMAX))
    BETA = 2# * (1# - Cos((JMAX - 1) * PI / JMAX))
    ALIM = Log(4# * JMAX / PI)
    K = 0
    Do
        K = K + 1
    Loop While 2 ^ K < ALIM
    EPS = 0.0001
    Call ADI(A(), B(), C(), D(), E(), F(), G(), U(), JMAX, K, ALPHA, BETA, EPS)
    Print
    Print Tab(5); "ADI Solution:"
    Print
    For I = 1 To JMAX
        For J = 1 To JMAX
            Print Tab(5 + (J - 1) * 8); Format$(U(I, J), "#.#0");
        Next J
    Next I
    Print Tab(5)
    Print Tab(5); "Test that sulotion satisfies Difference Eqns:"
    Print
    For I = 2 To JMAX - 1
        For J = 2 To JMAX - 1
            AAA = -4# * U(I, J) + U(I + 1, J)
            G(I, J) = AAA + U(I - 1, J) + U(I, J - 1) + U(I, J + 1)
        Next J
        For J = 2 To JMAX - 1
            Print Tab(5 + (J - 2) * 8); Format$(G(I, J), "#.#0");
        Next J
    Next I
End Sub
Sub ADI(A(), B(), C(), D(), E(), F(), G(), U(), JMAX, K, ALPHA, BETA, EPS)
    JJ = 100
    KK = 6
    NRR = 2 ^ (KK - 1)
    MAXITS = 100
    ZERO = 0#
    TWO = 2#
    HALF = 0.5
    Dim AA(100), BB(100), CC(100), RR(100), UU(100)
    Dim PSI(100, 100), ALPH(6), BET(6), R(32), S(32, 6)
    If JMAX > JJ Then Print " 'Increase JJ'"
    If K > KK - 1 Then Print " 'Increase KK'"
    K1 = K + 1
    NR = 2 ^ K
    ALPH(1) = ALPHA
    BET(1) = BETA
    For J = 1 To K
        ALPH(J + 1) = Sqr(ALPH(J) * BET(J))
        BET(J + 1) = HALF * (ALPH(J) + BET(J))
    Next J
    S(1, 1) = Sqr(ALPH(K1) * BET(K1))
    For J = 1 To K
        AB = ALPH(K1 - J) * BET(K1 - J)
        For N = 1 To 2 ^ (J - 1)
            DISC = Sqr(S(N, J) ^ 2 - AB)
            S(2 * N, J + 1) = S(N, J) + DISC
            S(2 * N - 1, J + 1) = AB / S(2 * N, J + 1)
        Next N
    Next J
    For N = 1 To NR
        R(N) = S(N, K1)
    Next N
    ANORMG = ZERO
    For J = 2 To JMAX - 1
      For L = 2 To JMAX - 1
        ANORMG = ANORMG + Abs(G(J, L))
        PSI(J, L) = -D(J, L) * U(J, L - 1) + (R(1) - E(J, L)) * U(J, L)
        PSI(J, L) = PSI(J, L) - F(J, L) * U(J, L + 1)
      Next L
    Next J
    NITS = MAXITS / NR
    For KITS = 1 To NITS
        For N = 1 To NR
            If N = NR Then
                NEXT1 = 1
            Else
                NEXT1 = N + 1
            End If
            RFACT = R(N) + R(NEXT1)
            For L = 2 To JMAX - 1
                For J = 2 To JMAX - 1
                    AA(J - 1) = A(J, L)
                    BB(J - 1) = B(J, L) + R(N)
                    CC(J - 1) = C(J, L)
                    RR(J - 1) = PSI(J, L) - G(J, L)
                Next J
                Call TRIDAG(AA(), BB(), CC(), RR(), UU(), JMAX - 2)
                For J = 2 To JMAX - 1
                    PSI(J, L) = -PSI(J, L) + TWO * R(N) * UU(J - 1)
                Next J
            Next L
            For J = 2 To JMAX - 1
                For L = 2 To JMAX - 1
                    AA(L - 1) = D(J, L)
                    BB(L - 1) = E(J, L) + R(N)
                    CC(L - 1) = F(J, L)
                    RR(L - 1) = PSI(J, L)
                Next L
                Call TRIDAG(AA, BB, CC, RR, UU, JMAX - 2)
                For L = 2 To JMAX - 1
                    U(J, L) = UU(L - 1)
                    PSI(J, L) = -PSI(J, L) + RFACT * UU(L - 1)
                Next L
            Next J
        Next N
        ANORM = ZERO
        For J = 2 To JMAX - 1
          For L = 2 To JMAX - 1
            RESID = A(J, L) * U(J - 1, L) + (B(J, L) + E(J, L)) * U(J, L)
            RESID = RESID + C(J, L) * U(J + 1, L)
            RESID = RESID + D(J, L) * U(J, L - 1) * U(J, L - 1)
            RESID = RESID + F(J, L) * U(J, L + 1) + G(J, L)
            ANROM = ANORM + Abs(RESID)
          Next L
        Next J
        If ANORM < EPS * ANORMG Then Exit Sub
    Next KITS
    Print " MAXITS exceeded"
End Sub
Sub TRIDAG(A(), B(), C(), R(), U(), N)
    NMAX = 100
    Dim GAM(100)
    If B(1) = 0# Then Exit Sub
    BET = B(1)
    U(1) = R(1) / BET
    For J = 2 To N
        GAM(J) = C(J - 1) / BET
        BET = B(J) - A(J) * GAM(J)
        If BET = 0# Then Exit Sub
        U(J) = (R(J) - A(J) * U(J - 1)) / BET
    Next J
    For J = N - 1 To 1 Step -1
        U(J) = U(J) - GAM(J + 1) * U(J + 1)
    Next J
End Sub