d15r1.txt

VB常用数值算法集2

Dim C2, FACTR, M, N
Private Sub Command1_Click()
    'PROGRAM D15R1
    'Driver for routine SHOOT
    'Solves for eigenvalues of Spheroidal Harmonics. Both
    'Prolate and Oblate case are handled simultaneously,
    'leading to six first-order equations.  Unknown to
    'SHOOT, these are actually two independent sets of
    'three coupled equations, one set with c^2 positive
    'and the other with c^2 negative.
    NVAR = 6: N2 = 2: DELTA = 0.001: EPS = 0.000001
    Dim V(2), DELV(2), F(2), DV(2)
    DX = 0.0001
    Do
      Print "     Input M,N,C-Squared (999 to end)"
      M = 2
      N = 2
      C2 = 0.1
      If C2 = 999 Then End
    Loop While (N < M) Or (M < 0) Or (N < 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
    V(1) = N * (N + 1) - M * (M + 1) + C2 / 2#
    V(2) = N * (N + 1) - M * (M + 1) - C2 / 2#
    DELV(1) = DELTA * V(1)
    DELV(2) = DELV(1)
    H1 = 0.1
    HMIN = 0#
    X1 = -1# + DX
    X2 = 0#
    Print "            Prolate                   Oblate"
    Print "     Mu(M,N)      Error Est.   Mu(M,N)     Error Est."
    Do
      Call SHOOT(NVAR, V(), DELV(), N2, X1, X2, EPS, H1, HMIN, F(), DV())
      Print Tab(5); Format$(V(1), "#.#####0"); Tab(18); Format$(DV(1), "#.#####0"),
      Print Tab(31); Format$(V(2), "#.#####0"), Tab(45); Format$(DV(2), "#.#####0")
    Loop While Abs(DV(1)) > Abs(EPS * V(1)) Or Abs(DV(2)) > Abs(EPS * V(2))
End Sub
Sub Load(X1, V(), Y())
    Y(3) = V(1)
    Y(2) = -(Y(3) - C2) * FACTR / 2# / (M + 1#)
    Y(1) = FACTR + Y(2) * DX
    Y(6) = V(2)
    Y(5) = -(Y(6) + C2) * FACTR / 2# / (M + 1#)
    Y(4) = FACTR + Y(5) * DX
End Sub
Sub SCORE(X2, Y(), F())
    If ((N - M) Mod 2) = 0 Then
        F(1) = Y(2)
        F(2) = Y(5)
    Else
        F(1) = Y(1)
        F(2) = Y(4)
    End If
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)
    DYDX(4) = Y(5)
    DYDX(6) = 0#
    DYDX(5) = (2# * X * (M + 1#) * Y(5) - (Y(6) + C2 * X * X) * Y(4)) / (1# - X * X)
End Sub