d15r2.frm

常用的数值算法的VB程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5325
   ClientLeft      =   2865
   ClientTop       =   1350
   ClientWidth     =   5820
   LinkTopic       =   "Form1"
   ScaleHeight     =   5325
   ScaleWidth      =   5820
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3840
      TabIndex        =   0
      Top             =   4560
      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 C2, M, N, FACTR, DX
Dim KMAX, KOUNT, DXSAV
Private Sub Command1_Click()
    'PROGRAM D15R2
    'Driver for routine SHOOTF
    NVAR = 3: N1 = 2: N2 = 1: DELTA = 0.001: EPS = 0.000001
    DXX = 0.0001: KMAX = 100
    Dim V1(1), DELV1(1), V2(2), DELV2(2), DV1(1), DV2(2), F(3)
    Do
      Print "Input M,N,C-Squared (999 to end)"
      M = 2
      N = 2
      C2 = 0.1
      If C2 = 999 Then End
      If (N < M) Or (M < 0) Then Print "Improper arguments"
    Loop While N < M Or M < 0
    Print Tab(5); M; Tab(10); N; Tab(15); C2
    FACTR = 1
    If M <> 0 Then
        Q1 = N
        For I = 1 To M
            FACTR = -0.5 * FACTR * (N + I) * (Q1 / I)
            Q1 = Q1 - 1
        Next I
    End If
    DX = DXX
    V1(1) = N * (N + 1) - M * (M + 1) + C2 / 2#
    If ((N - M) Mod 2) = 0 Then
        V2(1) = FACTR
    Else
        V2(1) = -FACTR
    End If
    V2(2) = V1(1) + 1
    DELV1(1) = DELTA * V1(1)
    DELV2(1) = DELTA * FACTR
    DELV2(2) = DELV1(1)
    H1 = 0.1
    HMIN = 0
    X1 = -1 + DX
    X2 = 1 - DX
    XF = 0
    Print "           Mu(-1)      Y(1-dx)          Mu(+1)"
    Do
    Call SHOOTF(NVAR, V1(), V2(), DELV1(), DELV2(), N1, N2, X1, X2, XF, EPS, H1, HMIN, F(), DV1(), DV2())
    Print Tab(2); "V       "; Format$(V1(1), "##.#####0"); Tab(22); Format$(V2(1), "#.#####0");
    Print Tab(40); Format$(V2(2), "##.#####0")
    Print Tab(2); "DV      "; Format$(DV1(1), "#.#####0"); Tab(22); Format$(DV2(1), "#.#####0");
    Print Tab(40); Format$(DV2(2), "#.#####0")
    Print
    Loop While Abs(DV1(1)) > Abs(EPS * V1(1))
End Sub
Sub LOAD1(X1, V1(), Y())
    Y(3) = V1(1)
    Y(2) = -(Y(3) - C2) * FACTR / 2# / (M + 1#)
    Y(1) = FACTR + Y(2) * DX
End Sub
Sub LOAD2(X2, V2(), Y())
    Y(3) = V2(2)
    Y(2) = (Y(3) - C2) * Y(1) / 2# / (M + 1#)
    Y(1) = V2(1)
End Sub
Sub SCORE(XF, Y(), F())
    For I = 1 To 3
        F(I) = Y(I)
    Next I
End Sub
Sub DERIVS(X, Y(), DYDX())
    DYDX(1) = Y(2)
    DYDX(3) = 0#
    DYDX(2) = (2# * X * (M + 1#) * Y(2) - (Y(3) - C2 * X * X) * Y(1)) / (1# - X * X)
End Sub
Sub SHOOTF(NVAR, V1(), V2(), DELV1(), DELV2(), N1, N2, X1, X2, XF, EPS, H1, HMIN, F(), DV1(), DV2())
    Dim Y(20), F1(20), F2(20), DFDV(20, 20), INDX(20), XP(200), YP(10, 200)
    Call LOAD1(X1, V1(), Y())
    Call ODEINT(Y(), NVAR, X1, XF, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    Call SCORE(XF, Y(), F1())
    Call LOAD2(X2, V2(), Y())
    Call ODEINT(Y(), NVAR, X2, XF, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    Call SCORE(XF, Y(), F2())
    J = 0
    For IV = 1 To N2
        J = J + 1
        SAV = V1(IV)
        V1(IV) = V1(IV) + DELV1(IV)
        Call LOAD1(X1, V1(), Y())
        Call ODEINT(Y(), NVAR, X1, XF, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
        Call SCORE(XF, Y(), F())
        For I = 1 To NVAR
            DFDV(I, J) = (F(I) - F1(I)) / DELV1(IV)
        Next I
        V1(IV) = SAV
    Next IV
    For IV = 1 To N1
        J = J + 1
        SAV = V2(IV)
        V2(IV) = V2(IV) + DELV2(IV)
        Call LOAD2(X2, V2(), Y())
        Call ODEINT(Y(), NVAR, X2, XF, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
        Call SCORE(XF, Y(), F())
        For I = 1 To NVAR
            DFDV(I, J) = (F2(I) - F(I)) / DELV2(IV)
        Next I
        V2(IV) = SAV
    Next IV
    For I = 1 To NVAR
        F(I) = F1(I) - F2(I)
        F1(I) = -F(I)
    Next I
    Call LUDCMP(DFDV(), NVAR, INDX(), DET)
    Call LUBKSB(DFDV(), NVAR, INDX(), F1())
    J = 0
    For IV = 1 To N2
        J = J + 1
        V1(IV) = V1(IV) + F1(J)
        DV1(IV) = F1(J)
    Next IV
    For IV = 1 To N1
        J = J + 1
        V2(IV) = V2(IV) + F1(J)
        DV2(IV) = F1(J)
    Next IV
    Erase INDX, DFDV, F2, F1, Y
End Sub
Sub ODEINT(YSTART(), NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    MAXSTP = 10000
    TWO = 2#
    ZERO = 0#
    TINY = 1E-30
    Dim YSCAL(10), Y(10), DYDX(10)
    KMAX = 100
    DXSAV = (X2 - X1) / 20
    X = X1
    H = Abs(H1) * Sgn(X2 - X1)
    NOK = 0
    NBAD = 0
    KOUNT = 0
    For I = 1 To NVAR
        Y(I) = YSTART(I)
    Next I
    If KMAX > 0 Then XSAV = X - DXSAV * TWO
    For NSTP = 1 To MAXSTP
        Call DERIVS(X, Y(), DYDX())
        For I = 1 To NVAR
            YSCAL(I) = Abs(Y(I)) + Abs(H * DYDX(I)) + TINY
        Next I
        If KMAX > 0 Then
            If Abs(X - XSAV) > Abs(DXSAV) Then
                If KOUNT < KMAX - 1 Then
                    KOUNT = KOUNT + 1
                    XP(KOUNT) = X
                    For I = 1 To NVAR
                        YP(I, KOUNT) = Y(I)
                    Next I
                    XSAV = X
                End If
            End If
        End If
        If (X + H - X2) * (X + H - X1) > ZERO Then H = X2 - X
        Call RKQC(Y(), DYDX(), NVAR, X, H, EPS, YSCAL(), HDID, HNEXT)
        If HDID = H Then
            NOK = NOK + 1
        Else
            NBAD = NBAD + 1
        End If
        If (X - X2) * (X2 - X1) >= ZERO Then
            For I = 1 To NVAR
                YSTART(I) = Y(I)
            Next I
            If KMAX <> 0 Then
                KOUNT = KOUNT + 1
                XP(KOUNT) = X
                For I = 1 To NVAR
                    YP(I, KOUNT) = Y(I)
                Next I
            End If
            Erase DYDX, Y, YSCAL
            Exit Sub
        End If
        If Abs(HNEXT) < HMIN Then
            Print " Stepsize smaller than minimum."
            Exit Sub
        End If
        H = HNEXT
    Next NSTP
    Print " Too many steps."
End Sub
Sub RKQC(Y(), DYDX(), N, X, HTRY, EPS, YSCAL(), HDID, HNEXT)
      ONE = 1#
      SAFETY = 0.9
      ERRCON = 0.0006
      FCOR = 0.066667
      Dim YTEMP(10), YSAV(10), DYSAV(10)
      PGROW = -0.2
      PSHRNK = -0.25
      XSAV = X
      For I = 1 To N
          YSAV(I) = Y(I)
          DYSAV(I) = DYDX(I)
      Next I
      H = HTRY
      Do
        HH = 0.5 * H
        Call RK4(YSAV(), DYSAV(), N, XSAV, HH, YTEMP())
        X = XSAV + HH
        Call DERIVS(X, YTEMP(), DYDX())
        Call RK4(YTEMP(), DYDX(), N, X, HH, Y())
        X = XSAV + H
        If X = XSAV Then
            Print " Stepsize not significant in RKQC."
            Exit Sub
        End If
        Call RK4(YSAV(), DYSAV(), N, XSAV, H, YTEMP())
        ERRMAX = 0#
        For I = 1 To N
            YTEMP(I) = Y(I) - YTEMP(I)
            If ERRMAX < Abs(YTEMP(I) / YSCAL(I)) Then
                ERRMAX = Abs(YTEMP(I) / YSCAL(I))
            End If
        Next I
        ERRMAX = ERRMAX / EPS
        If ERRMAX > ONE Then
            H = SAFETY * H * (ERRMAX ^ PSHRNK)
            FLAG = 1
        Else
            HDID = H
            If ERRMAX > ERRCON Then
                HNEXT = SAFETY * H * (ERRMAX ^ PGROW)
            Else
                HNEXT = 4# * H
            End If
            FLAG = 0
        End If
      Loop While FLAG = 1
      For I = 1 To N
          Y(I) = Y(I) + YTEMP(I) * FCOR
      Next I
      Erase DYSAV, YSAV, YTEMP
End Sub
Sub RK4(Y(), DYDX(), N, X, H, YOUT())
      Dim YT(10), DYT(10), DYM(10)
      HH = H * 0.5
      H6 = H / 6#
      XH = X + HH
      For I = 1 To N
          YT(I) = Y(I) + HH * DYDX(I)
      Next I
      Call DERIVS(XH, YT(), DYT())
      For I = 1 To N
          YT(I) = Y(I) + HH * DYT(I)
      Next I
      Call DERIVS(XH, YT(), DYM())
      For I = 1 To N
          YT(I) = Y(I) + H * DYM(I)
          DYM(I) = DYT(I) + DYM(I)
      Next I
      Call DERIVS(X + H, YT(), DYT())
      For I = 1 To N
          YOUT(I) = Y(I) + H6 * (DYDX(I) + DYT(I) + 2# * DYM(I))
      Next I
      Erase DYM, DYT, YT
End Sub
Sub LUBKSB(A(), N, INDX(), B())
    II = 0
    For I = 1 To N
        LL = INDX(I)
        Sum = B(LL)
        B(LL) = B(I)
        If II <> 0 Then
            For J = II To I - 1
                Sum = Sum - A(I, J) * B(J)
            Next J
        ElseIf Sum <> 0# Then
            II = I
        End If
        B(I) = Sum
    Next I
    For I = N To 1 Step -1
        Sum = B(I)
        If I < N Then
            For J = I + 1 To N
                Sum = Sum - A(I, J) * B(J)
            Next J
        End If
        B(I) = Sum / A(I, I)
    Next I
End Sub
Sub LUDCMP(A(), N, INDX(), D)
    NMAX = 100
    TINY = 1E-20
    Dim VV(100)
    D = 1#
    For I = 1 To N
        AAMAX = 0#
        For J = 1 To N
            If Abs(A(I, J)) > AAMAX Then AAMAX = Abs(A(I, J))
        Next J
        If AAMAX = 0# Then Print "Singular matrix."
        VV(I) = 1# / AAMAX
    Next I
    For J = 1 To N
        If J > 1 Then
            For I = 1 To J - 1
                Sum = A(I, J)
                If I > 1 Then
                    For K = 1 To I - 1
                        Sum = Sum - A(I, K) * A(K, J)
                    Next K
                    A(I, J) = Sum
                End If
            Next I
        End If
        AAMAX = 0#
        For I = J To N
            Sum = A(I, J)
            If J > 1 Then
                For K = 1 To J - 1
                    Sum = Sum - A(I, K) * A(K, J)
                Next K
                A(I, J) = Sum
            End If
            DUM = VV(I) * Abs(Sum)
            If DUM >= AAMAX Then
                IMAX = I
                AAMAX = DUM
            End If
        Next I
        If J <> IMAX Then
            For K = 1 To N
                DUM = A(IMAX, K)
                A(IMAX, K) = A(J, K)
                A(J, K) = DUM
            Next K
            D = -D
            VV(IMAX) = VV(J)
        End If
        INDX(J) = IMAX
        If J <> N Then
            If A(J, J) = 0# Then A(J, J) = TINY
            DUM = 1# / A(J, J)
            For I = J + 1 To N
                A(I, J) = A(I, J) * DUM
            Next I
        End If
    Next J
    If A(N, N) = 0# Then A(N, N) = TINY
End Sub