d12r4.frm

常用数值算法集哈

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5250
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   5250
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3000
      TabIndex        =   0
      Top             =   4440
      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 KMAX, KOUNT, DXSAV
Private Sub Command1_Click()
    'PROGRAM D12R4
    'Driver for routine ODEINT
    NVAR = 4
    Dim VSTART(4), XP(200), YP(10, 200)
    X1 = 1#
    X2 = 10#
    VSTART(1) = BESSJ0(X1)
    VSTART(2) = BESSJ1(X1)
    VSTART(3) = BESSJ(2, X1)
    VSTART(4) = BESSJ(3, X1)
    EPS = 0.0001
    H1 = 0.1
    HMIN = 0#
    KMAX = 100
    DXSAV = (X2 - X1) / 20#
    Call ODEINT(VSTART, NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    Print
    Print Tab(5); "Successful step: "; NOK
    Print Tab(5); "Bad step: "; NBAD
    Print Tab(5); "Stored intermediate values:  "; KOUNT
    Print
    Print Tab(5); "   X          Integral      BESSJ(3,X)"
    For I = 1 To KOUNT
        Print Tab(5); Format$(XP(I), "##.###0");
        Print Tab(19); Format$(YP(4, I), "##.#####0");
        Print Tab(33); Format$(BESSJ(3, XP(I)), "##.#####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 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)
    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
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