adi.txt

用VB语言实现偏微分方程的解法 包括解边值问题的松弛法、交替方向隐式方法等

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