📄 d12r6.frm
字号:
VERSION 5.00
Begin VB.Form Form1
Caption = "Form1"
ClientHeight = 4755
ClientLeft = 60
ClientTop = 345
ClientWidth = 5760
LinkTopic = "Form1"
ScaleHeight = 4755
ScaleWidth = 5760
StartUpPosition = 3 'Windows Default
Begin VB.CommandButton Command1
Caption = "Command1"
Height = 375
Left = 4080
TabIndex = 0
Top = 4080
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 D12R6
'Driver for routine BSSTEP
N = 4
Dim Y(4), DYDX(4), YSCAL(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)
For I = 1 To N
YSCAL(I) = 1#
Next I
HTRY = 0.1
Print
Print Tab(5); " eps htry hdid hnext"
For I = 1 To 15
EPS = Exp(-CSng(I))
Call BSSTEP(Y(), DYDX(), N, X, HTRY, EPS, YSCAL(), HDID, HNEXT)
Print Tab(5); Format$(EPS, ".###0E+00");
Print Tab(19); Format$(HTRY, ".##");
Print Tab(29); Format$(HDID, ".###000");
Print Tab(43); Format$(HNEXT, ".###000")
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 BSSTEP(Y(), DYDX(), NV, X1, HTRY, EPS, YSCAL(), HDID, HNEXT)
IMAX = 11
NUSE = 7
ONE = 1#
SHRINK = 0.95
GROW = 1.2
Dim YERR(10), YSAV(10), DYSAV(10), YSEQ(10), NSEQ(11)
NSEQ(1) = 2#: NSEQ(2) = 4#: NSEQ(3) = 6#: NSEQ(4) = 8#: NSEQ(5) = 12#
NSEQ(6) = 16#: NSEQ(7) = 24#: NSEQ(8) = 32#: NSEQ(9) = 48#
NSEQ(10) = 64#: NSEQ(11) = 96#
H = HTRY
XSAV = X1
For I = 1 To NV
YSAV(I) = Y(I)
DYSAV(I) = DYDX(I)
Next I
Do
For I = 1 To IMAX
Call MMID(YSAV(), DYSAV(), NV, XSAV, H, NSEQ(I), YSEQ())
XEST = (H / NSEQ(I)) ^ 2
Call RZEXTR(I, XEST, YSEQ(), Y(), YERR(), NV, NUSE)
If I > 3 Then
ERRMAX = 0#
For J = 1 To NV
If Abs(YERR(J) / YSCAL(J)) > ERRMAX Then
ERRMAX = Abs(YERR(J) / YSCAL(J))
End If
Next J
ERRMAX = ERRMAX / EPS
If ERRMAX < ONE Then
X1 = X1 + H
HDID = H
If I = NUSE Then
HNEXT = H * SHRINK
ElseIf I = NUSE - 1 Then
HNEXT = H * GROW
Else
HNEXT = (H * NSEQ(NUSE - 1)) / NSEQ(I)
End If
Erase NSEQ, YSEQ, DYSAV, YSAV, YERR
Exit Sub
End If
End If
Next I
H = 0.25 * H / 2 ^ ((IMAX - NUSE) / 2)
Loop Until X1 + H = X1
Print " Step size underflow"
End Sub
Sub RZEXTR(IEST, XEST, YEST(), YZ(), DY(), NV, NUSE)
Dim FX(7)
Static X(11), D(10, 7)
X(IEST) = XEST
If IEST = 1 Then
For J = 1 To NV
YZ(J) = YEST(J)
D(J, 1) = YEST(J)
DY(J) = YEST(J)
Next J
Else
M1 = IEST
If NUSE < IEST Then M1 = NUSE
For K = 1 To M1 - 1
FX(K + 1) = X(IEST - K) / XEST
Next K
For J = 1 To NV
YY = YEST(J)
V = D(J, 1)
C = YY
D(J, 1) = YY
For K = 2 To M1
B1 = FX(K) * V
B = B1 - C
If B <> 0# Then
B = (C - V) / B
DDY = C * B
C = B1 * B
Else
DDY = V
End If
If K <> M1 Then V = D(J, K)
D(J, K) = DDY
YY = YY + DDY
Next K
DY(J) = DDY
YZ(J) = YY
Next J
End If
Erase FX
End Sub
Sub MMID(Y(), DYDX(), NVAR, XS, HTOT, NSTEP, YOUT())
Dim YM(10), YN(10)
H = HTOT / NSTEP
For I = 1 To NVAR
YM(I) = Y(I)
YN(I) = Y(I) + H * DYDX(I)
Next I
X = XS + H
Call DERIVS(X, YN(), YOUT())
H2 = 2# * H
For N = 2 To NSTEP
For I = 1 To NVAR
SWAP = YM(I) + H2 * YOUT(I)
YM(I) = YN(I)
YN(I) = SWAP
Next I
X = X + H
Call DERIVS(X, YN(), YOUT())
Next N
For I = 1 To NVAR
YOUT(I) = 0.5 * (YM(I) + YN(I) + H * YOUT(I))
Next I
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -