?? sfroid.frm
字號(hào):
VERSION 5.00
Begin VB.Form Form1
Caption = "Form1"
ClientHeight = 8430
ClientLeft = 660
ClientTop = 345
ClientWidth = 7035
LinkTopic = "Form1"
ScaleHeight = 8430
ScaleWidth = 7035
Begin VB.CommandButton Command1
Caption = "Command1"
Height = 375
Left = 5640
TabIndex = 0
Top = 120
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 X(41), H, MM, N, C2, ANORM
Private Sub Command1_Click()
'PROGRAM SFROID
NE = 3: M = 41: NB = 1: NCI = 3: NCJ = 3: NCK = 42
NSI = 3: NSJ = 7: NYJ = 3: NYK = 41
Dim SCALV(3), INDEXV(3), Y(3, 41), C(3, 3, 42), S(3, 7)
ITMAX = 100
CONV = 0.000005
SLOWC = 1#
H = 1# / (M - 1)
C2 = 0#
Print "ENTER M,N"
MM = 2
N = 2
If (N + MM Mod 2) = 1 Then
INDEXV(1) = 1
INDEXV(2) = 2
INDEXV(3) = 3
Else
INDEXV(1) = 2
INDEXV(2) = 1
INDEXV(3) = 3
End If
ANORM = 1#
If MM <> 0 Then
Q1 = N
For I = 1 To MM
ANORM = -0.5 * ANORM * (N + I) * (Q1 / I)
Q1 = Q1 - 1#
Next I
End If
For K = 1 To M - 1
X(K) = (K - 1) * H
FAC1 = 1# - X(K) ^ 2
FAC2 = FAC1 ^ (-MM / 2#)
Y(1, K) = PLGNDR(N, MM, X(K)) * FAC2
DERIV = -((N - MM + 1) * PLGNDR(N + 1, MM, X(K)) - (N + 1) * X(K) * PLGNDR(N, MM, X(K))) / FAC1
Y(2, K) = MM * X(K) * Y(1, K) / FAC1 + DERIV * FAC2
Y(3, K) = N * (N + 1) - MM * (MM + 1)
Next K
X(M) = 1#
Y(1, M) = ANORM
Y(3, M) = N * (N + 1) - MM * (MM + 1)
Y(2, M) = (Y(3, M) - C2) * Y(1, M) / (2# * (MM + 1#))
SCALV(1) = Abs(ANORM)
If Y(2, M) > Abs(ANORM) Then
SCALV(2) = Y(2, M)
Else
SCALV(2) = Abs(ANORM)
End If
If Y(3, M) > 1 Then
SCALV(3) = Y(3, M)
Else
SCALV(3) = 1
End If
Do
Print "ENTER C^2 OR 999 TO END"
MSG1$ = "ENTER C^2 OR 999 TO END"
MSG$ = InputBox$(MSG1$, "INPUT C^2", "0.1")
C2 = Val(MSG$)
Print C2
If C2 = 999 Or MSG$ = "" Then Exit Do
Call SOLVDE(ITMAX, CONV, SLOWC, SCALV(), INDEXV(), NE, NB, M, Y(), NYJ, NYK, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Print "M = "; MM, " N = "; N; " C^2 = "; C2;
Print " LAMBDA = "; Y(3, 1) + MM * (MM + 1)
Print
Loop
End Sub
Sub SOLVDE(ITMAX, CONV, SLOWC, SCALV(), INDEXV(), NE, NB, M, Y(), NYJ, NYK, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Dim ERMAX(10), KMAX(10)
K1 = 1
K2 = M
NVARS = NE * M
J1 = 1
J2 = NB
J3 = NB + 1
J4 = NE
J5 = J4 + J1
J6 = J4 + J2
J7 = J4 + J3
J8 = J4 + J4
J9 = J8 + J1
IC1 = 1
IC2 = NE - NB
IC3 = IC2 + 1
IC4 = NE
JC1 = 1
JCF = IC3
For IT = 1 To ITMAX
K = K1
Call DIFEQ(K, K1, K2, J9, IC3, IC4, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
Call PINVS(IC3, IC4, J5, J9, JC1, K1, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
For K = K1 + 1 To K2
KP = K - 1
Call DIFEQ(K, K1, K2, J9, IC1, IC4, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
Call RED(IC1, IC4, J1, J2, J3, J4, J9, IC3, JC1, JCF, KP, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Call PINVS(IC1, IC4, J3, J9, JC1, K, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Next K
K = K2 + 1
Call DIFEQ(K, K1, K2, J9, IC1, IC2, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
Call RED(IC1, IC2, J5, J6, J7, J8, J9, IC3, JC1, JCF, K2, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Call PINVS(IC1, IC2, J7, J9, JCF, K2 + 1, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
Call BKSUB(NE, NB, JCF, K1, K2, C(), NCI, NCJ, NCK)
ERQ = 0#
For J = 1 To NE
JV = INDEXV(J)
ERMAX(J) = 0#
ERRJ = 0#
KMAX(J) = 0
VMAX = 0#
For K = K1 To K2
VZ = Abs(C(J, 1, K))
If VZ > VMAX Then
VMAX = VZ
KM = K
End If
ERRJ = ERRJ + VZ
Next K
ERQ = ERQ + ERRJ / SCALV(JV)
ERMAX(J) = C(J, 1, KM) / SCALV(JV)
KMAX(J) = KM
Next J
ERQ = ERQ / NVARS
If ERQ > SLOWC Then
DUM = ERQ
Else
DUM = SLOWC
End If
FAC = SLOWC / DUM
For JV = 1 To NE
J = INDEXV(JV)
For K = K1 To K2
Y(J, K) = Y(J, K) - FAC * C(JV, 1, K)
Next K
Next JV
Print Tab(3); Format$(IT, "####"); Tab(9); Format$(ERQ, "#.#####0");
Print Tab(19); Format$(FAC, "#.#####0")
Print
For J = 1 To NE
Print Tab(9); Format$(KMAX(J), "##0"); Tab(19); Format$(ERMAX(J), "#.#####0")
Next J
Print
If ERQ < CONV Then Exit For
Next IT
End Sub
Sub DIFEQ(K, K1, K2, JSF, IS1, ISF, INDEXV(), NE, S(), NSI, NSJ, Y(), NYJ, NYK)
M = 41
If K = K1 Then
If N + MM Mod 2 = 1 Then
S(3, 3 + INDEXV(1)) = 1# '方程(15-37)
S(3, 3 + INDEXV(2)) = 0#
S(3, 3 + INDEXV(3)) = 0#
S(3, JSF) = Y(1, 1) '方程(15-31)
Else
S(3, 3 + INDEXV(1)) = 0# '方程(15-37)
S(3, 3 + INDEXV(2)) = 1#
S(3, 3 + INDEXV(3)) = 0#
S(3, JSF) = Y(2, 1) '方程(15-31)
End If
ElseIf K > K2 Then
S(1, 3 + INDEXV(1)) = -(Y(3, M) - C2) / (2# * (MM + 1#)) '方程(15-38)
S(1, 3 + INDEXV(2)) = 1#
S(1, 3 + INDEXV(3)) = -Y(1, M) / (2# * (MM + 1#))
S(1, JSF) = Y(2, M) - (Y(3, M) - C2) * Y(1, M) / (2# * (MM + 1#)) '方程(15-32)
S(2, 3 + INDEXV(1)) = 1# '方程(15-39)
S(2, 3 + INDEXV(2)) = 0#
S(2, 3 + INDEXV(3)) = 0#
S(2, JSF) = Y(1, M) - ANORM '方程(15-33)
Else
S(1, INDEXV(1)) = -1# '方程(15-34)
S(1, INDEXV(2)) = -0.5 * H
S(1, INDEXV(3)) = 0#
S(1, 3 + INDEXV(1)) = 1#
S(1, 3 + INDEXV(2)) = -0.5 * H
S(1, 3 + INDEXV(3)) = 0#
TEMP = H / (1# - (X(K) + X(K - 1)) ^ 2 * 0.25)
TEMP2 = 0.5 * (Y(3, K) + Y(3, K - 1)) - C2 * 0.25 * (X(K) + X(K - 1)) ^ 2
S(2, INDEXV(1)) = TEMP * TEMP2 * 0.5 '方程(15-35)
S(2, INDEXV(2)) = -1# - 0.5 * TEMP * (MM + 1#) * (X(K) + X(K - 1))
S(2, INDEXV(3)) = 0.25 * TEMP * (Y(1, K) + Y(1, K - 1))
S(2, 3 + INDEXV(1)) = S(2, INDEXV(1))
S(2, 3 + INDEXV(2)) = 2# + S(2, INDEXV(2))
S(2, 3 + INDEXV(3)) = S(2, INDEXV(3))
S(3, INDEXV(1)) = 0# '方程(15-36)
S(3, INDEXV(2)) = 0#
S(3, INDEXV(3)) = -1#
S(3, 3 + INDEXV(1)) = 0#
S(3, 3 + INDEXV(2)) = 0#
S(3, 3 + INDEXV(3)) = 1#
S(1, JSF) = Y(1, K) - Y(1, K - 1) - 0.5 * H * (Y(2, K) + Y(2, K - 1)) '方程(15-26)
DUM = (X(K) + X(K - 1)) * 0.5 * (MM + 1#) * (Y(2, K) + Y(2, K - 1)) '方程(15-27)
DUM = DUM - TEMP2 * 0.5 * (Y(1, K) + Y(1, K - 1))
S(2, JSF) = Y(2, K) - Y(2, K - 1) - TEMP * DUM
S(3, JSF) = Y(3, K) - Y(3, K - 1) '方程(15-30)
End If
End Sub
Sub BKSUB(NE, NB, JF, K1, K2, C(), NCI, NCJ, NCK)
NBF = NE - NB
For K = K2 To K1 Step -1
KP = K + 1
For J = 1 To NBF
XX = C(J, JF, KP)
For I = 1 To NE
C(I, JF, K) = C(I, JF, K) - C(I, J, K) * XX
Next I
Next J
Next K
For K = K1 To K2
KP = K + 1
For I = 1 To NB
C(I, 1, K) = C(I + NBF, JF, K)
Next I
For I = 1 To NBF
C(I + NB, 1, K) = C(I, JF, KP)
Next I
Next K
End Sub
Sub PINVS(IE1, IE2, JE1, JSF, JC1, K, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
ZERO = 0#: ONE = 1
Dim PSCL(10), INDXR(10)
JE2 = JE1 + IE2 - IE1
JS1 = JE2 + 1
For I = IE1 To IE2
BIG = ZERO
For J = JE1 To JE2
If Abs(S(I, J)) > BIG Then BIG = Abs(S(I, J))
Next J
If BIG = ZERO Then Print "Singular matrix, row all 0"
PSCL(I) = ONE / BIG
INDXR(I) = 0
Next I
For ID = IE1 To IE2
PIV = ZERO
For I = IE1 To IE2
If INDXR(I) = 0 Then
BIG = ZERO
For J = JE1 To JE2
If Abs(S(I, J)) > BIG Then
JP = J
BIG = Abs(S(I, J))
End If
Next J
If BIG * PSCL(I) > PIV Then
IPIV = I
JPIV = JP
PIV = BIG * PSCL(I)
End If
End If
Next I
If S(IPIV, JPIV) = ZERO Then Print "Singular matrix"
INDXR(IPIV) = JPIV
PIVINV = ONE / S(IPIV, JPIV)
For J = JE1 To JSF
S(IPIV, J) = S(IPIV, J) * PIVINV
Next J
S(IPIV, JPIV) = ONE
For I = IE1 To IE2
If INDXR(I) <> JPIV Then
If S(I, JPIV) <> ZERO Then
DUM = S(I, JPIV)
For J = JE1 To JSF
S(I, J) = S(I, J) - DUM * S(IPIV, J)
Next J
S(I, JPIV) = ZERO
End If
End If
Next I
Next ID
JCOFF = JC1 - JS1
ICOFF = IE1 - JE1
For I = IE1 To IE2
IROW = INDXR(I) + ICOFF
For J = JS1 To JSF
C(IROW, J + JCOFF, K) = S(I, J)
Next J
Next I
End Sub
Sub RED(IZ1, IZ2, JZ1, JZ2, JM1, JM2, JMF, IC1, JC1, JCF, KC, C(), NCI, NCJ, NCK, S(), NSI, NSJ)
LOFF = JC1 - JM1
IC = IC1
For J = JZ1 To JZ2
For L = JM1 To JM2
VX = C(IC, L + LOFF, KC)
For I = IZ1 To IZ2
S(I, L) = S(I, L) - S(I, J) * VX
Next I
Next L
VX = C(IC, JCF, KC)
For I = IZ1 To IZ2
S(I, JMF) = S(I, JMF) - S(I, J) * VX
Next I
IC = IC + 1
Next J
End Sub
Function PLGNDR(L, M, X)
If M < 0 Or M > L Or Abs(X) > 1# Then Print "bad arguments"
PMM = 1#
If M > 0 Then
SOMX2 = Sqr((1# - X) * (1# + X))
FACT = 1#
For I = 1 To M
PMM = -PMM * FACT * SOMX2
FACT = FACT + 2#
Next I
End If
If L = M Then
PLGNDR = PMM
Else
PMMP1 = X * (2 * M + 1) * PMM
If L = M + 1 Then
PLGNDR = PMMP1
Else
For LL = M + 2 To L
PLL = X * (2 * LL - 1) * PMMP1 - (LL + M - 1) * PMM
PLL = PLL / (LL - M)
PMM = PMMP1
PMMP1 = PLL
Next LL
PLGNDR = PLL
End If
End If
End Function
?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -