d12r1.frm

常用数值算法集哈

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5655
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   6030
   LinkTopic       =   "Form1"
   ScaleHeight     =   5655
   ScaleWidth      =   6030
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   4320
      TabIndex        =   0
      Top             =   5040
      Width           =   1215
   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 D12R1
    'Driver for routine RK4
    N = 4
    Dim Y(4), DYDX(4), YOUT(4)
    X = 1#
    Y(1) = BESSJ0(X)
    Y(2) = BESSJ1(X)
    Y(3) = BESSJ(2, X)
    Y(4) = BESSJ(3, X)
    DYDX(1) = -Y(2)
    DYDX(2) = Y(1) - Y(2)
    DYDX(3) = Y(2) - 2# * Y(3)
    DYDX(4) = Y(3) - 3# * Y(4)
    Print Tab(5); "Bessel Function:   J0   J1   J2   J3"
    For I = 1 To 5
        H = 0.2 * I
        Call RK4(Y(), DYDX(), N, X, H, YOUT())
        Print
        Print Tab(5); "For a step size of: "; Format$(H, "#.#0")
        Print Tab(5); "RK4:      ";
        For J = 1 To 4
            Print Format$(YOUT(J), ".#####0"),
        Next J
        Print Tab(5); "Actual:",
        Print Tab(15); Format$(BESSJ0(X + H), ".#####0");
        Print Tab(29); Format$(BESSJ1(X + H), ".#####0");
        Print Tab(43); Format$(BESSJ(2, X + H), ".#####0");
        Print Tab(57); Format$(BESSJ(3, X + H), ".#####0")
    Next I
End Sub
Sub DERIVS(X, Y(), DYDX())
    DYDX(1) = -Y(2)
    DYDX(2) = Y(1) - (1# / X) * Y(2)
    DYDX(3) = Y(2) - (2# / X) * Y(3)
    DYDX(4) = Y(3) - (3# / X) * Y(4)
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
Function BESSJ(N, X)
    IACC = 40
    BIGNO = 10000000000#
    BIGNI = 0.0000000001
    If N < 2 Then
        Form1.Print "bad argument N in BASSJ"
        Exit Function
    End If
    AX = Abs(X)
    If AX = 0 Then
        BESSJ = 0#
    ElseIf AX > CSng(N) Then
        TOX = 2# / AX
        BJM = BESSJ0(AX)
        BJ = BESSJ1(AX)
        For J = 1 To N - 1
            BJP = J * TOX * BJ - BJM
            BJM = BJ
            BJ = BJP
        Next J
        BESSJ = BJ
    Else
        TOX = 2# / AX
        M = 2 * Int(((N + Int(Sqr(IACC * N)))) / 2)
        BESJ = 0#
        JSUM = 0
        Sum = 0#
        BJP = 0#
        BJ = 1#
        For J = M To 1 Step -1
            BJM = J * TOX * BJ - BJP
            BJP = BJ
            BJ = BJM
            If Abs(BJ) > BIGNO Then
                BJ = BJ * BIGNI
                BJP = BJP * BIGNI
                BESJ = BESJ * BIGNI
                Sum = Sum * BIGNI
            End If
            If JSUM <> 0 Then Sum = Sum + BJ
            JSUM = 1 - JSUM
            If J = N Then BESJ = BJP
        Next J
        Sum = 2# * Sum - BJ
        BESSJ = BESJ / Sum
    End If
End Function
Function BESSJ0(X)
    P1# = 1#
    P2# = -0.001098628627
    P3# = 0.00002734510407
    P4# = -0.000002073370639
    P5# = 2.093887211E-07
    Q1# = -0.01562499995
    Q2# = 0.0001430488765
    Q3# = -0.000006911147651
    Q4# = 7.621095161E-07
    Q5# = -9.34945152E-08
    R1# = 57568490574#
    R2# = -13362590354#
    R3# = 651619640.7
    R4# = -11214424.18
    R5# = 77392.33017
    R6# = -184.9052456
    S1# = 57568490411#
    S2# = 1029532985#
    S3# = 9494680.718
    S4# = 59272.64853
    S5# = 267.8532712
    S6# = 1#
    If Abs(X) < 8# Then
       Y# = X * X
       BBB# = Y# * (R4# + Y# * (R5# + Y# * R6#))
       AAA# = R1# + Y# * (R2# + Y# * (R3# + BBB#))
       CCC# = Y# * (S3# + Y# * (S4# + Y# * (S5# + Y# * S6#)))
       BESSJ0 = AAA / (S1# + Y# * (S2# + CCC#))
    Else
       AX# = Abs(X)
       Z# = 8# / AX#
       Y# = Z# * Z#
       XX = AX# - 0.785398164
       CCC# = Y# * (P3# + Y# * (P4# + Y# * P5#))
       AAA# = P1# + Y# * (P2# + CCC#)
       DDD# = Y# * (Q3# + Y# * (Q4# + Y# * Q5#))
       EEE# = Z# * Sin(XX) * (Q1# + Y# * (Q2# + DDD#))
       BESSJ0 = Sqr(0.636619772 / AX#) * (Cos(XX) * AAA# - EEE#)
    End If
End Function
Function BESSJ1(X)
    R1 = 72362614232#
    R2 = -7895059235#
    R3 = 242396853.1
    R4 = -2972611.439
    R5 = 15704.4826
    R6 = -30.16036606
    S1 = 144725228442#
    S2 = 2300535178#
    S3 = 18583304.74
    S4 = 99447.43394
    S5 = 376.9991397
    S6 = 1#
    P1 = 1#
    P2 = 0.00183105
    P3 = -0.00003516396496
    P4 = 0.000002457520174
    P5 = -0.000000240337019
    Q1 = 0.04687499995
    Q2 = -0.0002002690873
    Q3 = 0.000008449199096
    Q4 = -0.00000088228987
    Q5 = 0.000000105787412
    If Abs(X) < 8# Then
        Y# = X ^ 2
        AAA# = R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))))
        BBB# = S1 + Y * (S2 + Y * (S3 + Y * (S4 + Y * (S5 + Y * S6))))
        BESSJ1 = X * AAA# / BBB#
    Else
        AX = Abs(X)
        Z = 8# / AX
        Y# = Z ^ 2
        XX = AX - 2.356194491
        AAA# = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
        BBB# = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
        CCC# = Sqr(0.636619772 / AX)
        BESSJ1 = CCC# * (Cos(XX) * AAA# - Z * Sin(XX) * BBB# * Sgn(X))
    End If
End Function