sfroid.frm

常用的数值算法的VB程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   8430
   ClientLeft      =   660
   ClientTop       =   345
   ClientWidth     =   7035
   LinkTopic       =   "Form1"
   ScaleHeight     =   8430
   ScaleWidth      =   7035
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   5640
      TabIndex        =   0
      Top             =   120
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Dim X(41), H, MM, N, C2, ANORM
Private Sub Command1_Click()
    'PROGRAM SFROID
    NE = 3: M = 41: NB = 1: NCI = 3: NCJ = 3: NCK = 42
    NSI = 3: NSJ = 7: NYJ = 3: NYK = 41
    Dim SCALV(3), INDEXV(3), Y(3, 41), C(3, 3, 42), S(3, 7)
    ITMAX = 100
    CONV = 0.000005
    SLOWC = 1#
    H = 1# / (M - 1)
    C2 = 0#
    Print "ENTER M,N"
    MM = 2
    N = 2
    If (N + MM Mod 2) = 1 Then
        INDEXV(1) = 1
        INDEXV(2) = 2
        INDEXV(3) = 3
    Else
        INDEXV(1) = 2
        INDEXV(2) = 1
        INDEXV(3) = 3
    End If
    ANORM = 1#
    If MM <> 0 Then
        Q1 = N
        For I = 1 To MM
            ANORM = -0.5 * ANORM * (N + I) * (Q1 / I)
            Q1 = Q1 - 1#
        Next I
    End If
    For K = 1 To M - 1
        X(K) = (K - 1) * H
        FAC1 = 1# - X(K) ^ 2
        FAC2 = FAC1 ^ (-MM / 2#)
        Y(1, K) = PLGNDR(N, MM, X(K)) * FAC2
        DERIV = -((N - MM + 1) * PLGNDR(N + 1, MM, X(K)) - (N + 1) * X(K) * PLGNDR(N, MM, X(K))) / FAC1
        Y(2, K) = MM * X(K) * Y(1, K) / FAC1 + DERIV * FAC2
        Y(3, K) = N * (N + 1) - MM * (MM + 1)
    Next K
    X(M) = 1#
    Y(1, M) = ANORM
    Y(3, M) = N * (N + 1) - MM * (MM + 1)
    Y(2, M) = (Y(3, M) - C2) * Y(1, M) / (2# * (MM + 1#))
    SCALV(1) = Abs(ANORM)
    If Y(2, M) > Abs(ANORM) Then
        SCALV(2) = Y(2, M)
    Else
        SCALV(2) = Abs(ANORM)
    End If
    If Y(3, M) > 1 Then
        SCALV(3) = Y(3, M)
    Else
        SCALV(3) = 1
    End If
    Do
      Print "ENTER C^2 OR  999 TO END"
      MSG1$ = "ENTER C^2 OR  999 TO END"
      MSG$ = InputBox$(MSG1$, "INPUT C^2", "0.1")
      C2 = Val(MSG$)
      Print C2
      If C2 = 999 Or MSG$ = "" Then Exit Do
      Call SOLVDE(ITMAX, CONV, SLOWC, SCALV(), INDEXV(), NE, NB, M, Y(), NYJ, NYK, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
      Print "M = "; MM, "    N = "; N; "    C^2 = "; C2;
      Print "    LAMBDA = "; Y(3, 1) + MM * (MM + 1)
      Print
    Loop
End Sub
Sub SOLVDE(ITMAX, CONV, SLOWC, SCALV(), INDEXV(), NE, NB, M, Y(), NYJ, NYK, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
    Dim ERMAX(10), KMAX(10)
    K1 = 1
    K2 = M
    NVARS = NE * M
    J1 = 1
    J2 = NB
    J3 = NB + 1
    J4 = NE
    J5 = J4 + J1
    J6 = J4 + J2
    J7 = J4 + J3
    J8 = J4 + J4
    J9 = J8 + J1
    IC1 = 1
    IC2 = NE - NB
    IC3 = IC2 + 1
    IC4 = NE
    JC1 = 1
    JCF = IC3
    For IT = 1 To ITMAX
        K = K1
        Call DIFEQ(K, K1, K2, J9, IC3, IC4, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
        Call PINVS(IC3, IC4, J5, J9, JC1, K1, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
        For K = K1 + 1 To K2
            KP = K - 1
            Call DIFEQ(K, K1, K2, J9, IC1, IC4, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
            Call RED(IC1, IC4, J1, J2, J3, J4, J9, IC3, JC1, JCF, KP, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
            Call PINVS(IC1, IC4, J3, J9, JC1, K, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
        Next K
        K = K2 + 1
        Call DIFEQ(K, K1, K2, J9, IC1, IC2, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
        Call RED(IC1, IC2, J5, J6, J7, J8, J9, IC3, JC1, JCF, K2, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
        Call PINVS(IC1, IC2, J7, J9, JCF, K2 + 1, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
        Call BKSUB(NE, NB, JCF, K1, K2, C(), NCI, NCJ, NCK)
        ERQ = 0#
        For J = 1 To NE
            JV = INDEXV(J)
            ERMAX(J) = 0#
            ERRJ = 0#
            KMAX(J) = 0
            VMAX = 0#
            For K = K1 To K2
                VZ = Abs(C(J, 1, K))
                If VZ > VMAX Then
                    VMAX = VZ
                    KM = K
                End If
                ERRJ = ERRJ + VZ
            Next K
            ERQ = ERQ + ERRJ / SCALV(JV)
            ERMAX(J) = C(J, 1, KM) / SCALV(JV)
            KMAX(J) = KM
        Next J
        ERQ = ERQ / NVARS
        If ERQ > SLOWC Then
            DUM = ERQ
        Else
            DUM = SLOWC
        End If
        FAC = SLOWC / DUM
        For JV = 1 To NE
            J = INDEXV(JV)
            For K = K1 To K2
                Y(J, K) = Y(J, K) - FAC * C(JV, 1, K)
            Next K
        Next JV
        Print Tab(3); Format$(IT, "####"); Tab(9); Format$(ERQ, "#.#####0");
        Print Tab(19); Format$(FAC, "#.#####0")
        Print
        For J = 1 To NE
            Print Tab(9); Format$(KMAX(J), "##0"); Tab(19); Format$(ERMAX(J), "#.#####0")
        Next J
        Print
        If ERQ < CONV Then Exit For
    Next IT
End Sub
Sub DIFEQ(K, K1, K2, JSF, IS1, ISF, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
    M = 41
    If K = K1 Then
        If N + MM Mod 2 = 1 Then
            S(3, 3 + INDEXV(1)) = 1#  '方程(15-37)
            S(3, 3 + INDEXV(2)) = 0#
            S(3, 3 + INDEXV(3)) = 0#
            S(3, JSF) = Y(1, 1)       '方程(15-31)
        Else
            S(3, 3 + INDEXV(1)) = 0#  '方程(15-37)
            S(3, 3 + INDEXV(2)) = 1#
            S(3, 3 + INDEXV(3)) = 0#
            S(3, JSF) = Y(2, 1)       '方程(15-31)
        End If
    ElseIf K > K2 Then
        S(1, 3 + INDEXV(1)) = -(Y(3, M) - C2) / (2# * (MM + 1#))  '方程(15-38)
        S(1, 3 + INDEXV(2)) = 1#
        S(1, 3 + INDEXV(3)) = -Y(1, M) / (2# * (MM + 1#))
        S(1, JSF) = Y(2, M) - (Y(3, M) - C2) * Y(1, M) / (2# * (MM + 1#)) '方程(15-32)
        S(2, 3 + INDEXV(1)) = 1#  '方程(15-39)
        S(2, 3 + INDEXV(2)) = 0#
        S(2, 3 + INDEXV(3)) = 0#
        S(2, JSF) = Y(1, M) - ANORM  '方程(15-33)
    Else
        S(1, INDEXV(1)) = -1#  '方程(15-34)
        S(1, INDEXV(2)) = -0.5 * H
        S(1, INDEXV(3)) = 0#
        S(1, 3 + INDEXV(1)) = 1#
        S(1, 3 + INDEXV(2)) = -0.5 * H
        S(1, 3 + INDEXV(3)) = 0#
        TEMP = H / (1# - (X(K) + X(K - 1)) ^ 2 * 0.25)
        TEMP2 = 0.5 * (Y(3, K) + Y(3, K - 1)) - C2 * 0.25 * (X(K) + X(K - 1)) ^ 2
        S(2, INDEXV(1)) = TEMP * TEMP2 * 0.5  '方程(15-35)
        S(2, INDEXV(2)) = -1# - 0.5 * TEMP * (MM + 1#) * (X(K) + X(K - 1))
        S(2, INDEXV(3)) = 0.25 * TEMP * (Y(1, K) + Y(1, K - 1))
        S(2, 3 + INDEXV(1)) = S(2, INDEXV(1))
        S(2, 3 + INDEXV(2)) = 2# + S(2, INDEXV(2))
        S(2, 3 + INDEXV(3)) = S(2, INDEXV(3))
        S(3, INDEXV(1)) = 0#  '方程(15-36)
        S(3, INDEXV(2)) = 0#
        S(3, INDEXV(3)) = -1#
        S(3, 3 + INDEXV(1)) = 0#
        S(3, 3 + INDEXV(2)) = 0#
        S(3, 3 + INDEXV(3)) = 1#
        S(1, JSF) = Y(1, K) - Y(1, K - 1) - 0.5 * H * (Y(2, K) + Y(2, K - 1))  '方程(15-26)
        DUM = (X(K) + X(K - 1)) * 0.5 * (MM + 1#) * (Y(2, K) + Y(2, K - 1))    '方程(15-27)
        DUM = DUM - TEMP2 * 0.5 * (Y(1, K) + Y(1, K - 1))
        S(2, JSF) = Y(2, K) - Y(2, K - 1) - TEMP * DUM
        S(3, JSF) = Y(3, K) - Y(3, K - 1)  '方程(15-30)
    End If
End Sub
Sub BKSUB(NE, NB, JF, K1, K2, C(), NCI, NCJ, NCK)
    NBF = NE - NB
    For K = K2 To K1 Step -1
        KP = K + 1
        For J = 1 To NBF
            XX = C(J, JF, KP)
            For I = 1 To NE
                C(I, JF, K) = C(I, JF, K) - C(I, J, K) * XX
            Next I
        Next J
    Next K
    For K = K1 To K2
        KP = K + 1
        For I = 1 To NB
            C(I, 1, K) = C(I + NBF, JF, K)
        Next I
        For I = 1 To NBF
            C(I + NB, 1, K) = C(I, JF, KP)
        Next I
    Next K
End Sub
Sub PINVS(IE1, IE2, JE1, JSF, JC1, K, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
    ZERO = 0#: ONE = 1
    Dim PSCL(10), INDXR(10)
    JE2 = JE1 + IE2 - IE1
    JS1 = JE2 + 1
    For I = IE1 To IE2
        BIG = ZERO
        For J = JE1 To JE2
            If Abs(S(I, J)) > BIG Then BIG = Abs(S(I, J))
        Next J
        If BIG = ZERO Then Print "Singular matrix, row all 0"
        PSCL(I) = ONE / BIG
        INDXR(I) = 0
    Next I
    For ID = IE1 To IE2
        PIV = ZERO
        For I = IE1 To IE2
            If INDXR(I) = 0 Then
                BIG = ZERO
                For J = JE1 To JE2
                    If Abs(S(I, J)) > BIG Then
                        JP = J
                        BIG = Abs(S(I, J))
                    End If
                Next J
                If BIG * PSCL(I) > PIV Then
                    IPIV = I
                    JPIV = JP
                    PIV = BIG * PSCL(I)
                End If
            End If
        Next I
        If S(IPIV, JPIV) = ZERO Then Print "Singular matrix"
        INDXR(IPIV) = JPIV
        PIVINV = ONE / S(IPIV, JPIV)
        For J = JE1 To JSF
            S(IPIV, J) = S(IPIV, J) * PIVINV
        Next J
        S(IPIV, JPIV) = ONE
        For I = IE1 To IE2
            If INDXR(I) <> JPIV Then
                If S(I, JPIV) <> ZERO Then
                    DUM = S(I, JPIV)
                    For J = JE1 To JSF
                        S(I, J) = S(I, J) - DUM * S(IPIV, J)
                    Next J
                    S(I, JPIV) = ZERO
                End If
            End If
        Next I
    Next ID
    JCOFF = JC1 - JS1
    ICOFF = IE1 - JE1
    For I = IE1 To IE2
        IROW = INDXR(I) + ICOFF
        For J = JS1 To JSF
            C(IROW, J + JCOFF, K) = S(I, J)
        Next J
    Next I
End Sub
Sub RED(IZ1, IZ2, JZ1, JZ2, JM1, JM2, JMF, IC1, JC1, JCF, KC, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
    LOFF = JC1 - JM1
    IC = IC1
    For J = JZ1 To JZ2
        For L = JM1 To JM2
            VX = C(IC, L + LOFF, KC)
            For I = IZ1 To IZ2
                S(I, L) = S(I, L) - S(I, J) * VX
            Next I
        Next L
        VX = C(IC, JCF, KC)
        For I = IZ1 To IZ2
            S(I, JMF) = S(I, JMF) - S(I, J) * VX
        Next I
        IC = IC + 1
    Next J
End Sub
Function PLGNDR(L, M, X)
    If M < 0 Or M > L Or Abs(X) > 1# Then Print "bad arguments"
    PMM = 1#
    If M > 0 Then
        SOMX2 = Sqr((1# - X) * (1# + X))
        FACT = 1#
        For I = 1 To M
            PMM = -PMM * FACT * SOMX2
            FACT = FACT + 2#
        Next I
    End If
    If L = M Then
        PLGNDR = PMM
    Else
        PMMP1 = X * (2 * M + 1) * PMM
        If L = M + 1 Then
            PLGNDR = PMMP1
        Else
            For LL = M + 2 To L
                PLL = X * (2 * LL - 1) * PMMP1 - (LL + M - 1) * PMM
                PLL = PLL / (LL - M)
                PMM = PMMP1
                PMMP1 = PLL
            Next LL
            PLGNDR = PLL
        End If
    End If
End Function