d10r10.frm

VB常用数值算法集 内含有解线性代数方程组

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3060
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   3060
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      TabIndex        =   0
      Top             =   2400
      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 D10R10
    'Driver for routine LAGUER
    M = 4
    MP1 = M + 1
    NTRY = 21
    EPS = 0.000001
    Dim A(2, 5), Y(2, 21), X(2)
    For J = 1 To MP1
        For I = 1 To 2
            A(I, J) = 0
        Next I
    Next J
    A(2, 1) = 2#
    A(1, 3) = -1#
    A(2, 3) = -2#
    A(1, 5) = 1#
    Print
    Print Tab(5); "Roots of polynomial x^4-(1+2i)*x^2+2i"
    Print
    Print Tab(5); "             Real         Complex"
    Print
    N = 0
    POLISH% = 0
    For I = 1 To NTRY
        X(1) = (I - 11#) / 10#
        X(2) = (I - 11#) / 10#
        Call LAGUER(A(), M, X(), EPS, POLISH%)
        If N = 0 Then
            N = 1
            Y(1, 1) = X(1)
            Y(2, 1) = X(2)
            Print Tab(5); Format$(N, "#");
            Print Tab(15); Format$(X(1), "#.#####0");
            Print Tab(29); Format$(X(2), "#.#####0")
        Else
            IFLAG = 0
            For J = 1 To N
                AAA = CABS(X(1) - Y(1, J), X(2) - Y(2, J))
                BBB = EPS * CABS(X(1), X(2))
                If AAA <= BBB Then IFLAG = 1
            Next J
            If IFLAG = 0 Then
                N = N + 1
                Y(1, N) = X(1)
                Y(2, N) = X(2)
                Print Tab(5); Format$(N, "#");
                Print Tab(15); Format$(X(1), "#.#####0");
                Print Tab(29); Format$(X(2), "#.#####0")
            End If
        End If
    Next I
End Sub
Sub LAGUER(A(), M, X(), EPS, POLISH%)
    Dim ZERO(2), B(2), D(2), F(2), G(2), H(2)
    Dim G2(2), SQ(2), GP(2), GM(2), DX(2), X1(2)
    ZERO(1) = 0#
    ZERO(2) = 0#
    EPSS = 0.00000006
    MAXIT = 100
    DXOLD = CABS(X(1), X(2))
    For ITER = 1 To MAXIT
        B(1) = A(1, M + 1)
        B(2) = A(2, M + 1)
        ERQ = CABS(B(1), X(2))
        D(1) = ZERO(1)
        D(2) = ZERO(2)
        F(1) = ZERO(1)
        F(2) = ZERO(2)
        ABX = CABS(X(1), X(2))
        For J = M To 1 Step -1
            DUM = X(1) * F(1) - X(2) * F(2) + D(1)
            F(2) = X(2) * F(1) + X(1) * F(2) + D(2)
            F(1) = DUM
            DUM = X(1) * D(1) - X(2) * D(2) + B(1)
            D(2) = X(2) * D(1) + X(1) * D(2) + B(2)
            D(1) = DUM
            DUM = X(1) * B(1) - X(2) * B(2) + A(1, J)
            B(2) = X(2) * B(1) + X(1) * B(2) + A(2, J)
            B(1) = DUM
            ERQ = CABS(B(1), B(2)) + ABX * ERQ
        Next J
        ERQ = EPSS * ERQ
        If CABS(B(1), B(2)) <= ERQ Then
            Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
            Exit Sub
        Else
            G(1) = CDIV1(D(1), D(2), B(1), B(2))
            G(2) = CDIV2(D(1), D(2), B(1), B(2))
            G2(1) = G(1) * G(1) - G(2) * G(2)
            G2(2) = 2 * G(1) * G(2)
            H(1) = G2(1) - 2 * CDIV1(F(1), F(2), B(1), B(2))
            H(2) = G2(2) - 2 * CDIV2(F(1), F(2), B(1), B(2))
            DUM1 = (M - 1) * (M * H(1) - G2(1))
            DUM2 = (M - 1) * (M * H(2) - G2(2))
            SQ(1) = CSQR1(DUM1, DUM2)
            SQ(2) = CSQR2(DUM1, DUM2)
            GP(1) = G(1) + SQ(1)
            GP(2) = G(2) + SQ(2)
            GM(1) = G(1) - SQ(1)
            GM(2) = G(2) - SQ(2)
            If CABS(GP(1), GP(2)) < CABS(GM(1), GM(2)) Then
                GP(1) = GM(1)
                GP(2) = GM(2)
            End If
            DX(1) = CDIV1(M, 0, GP(1), GP(2))
            DX(2) = CDIV2(M, 0, GP(1), GP(2))
        End If
        X1(1) = X(1) - DX(1)
        X1(2) = X(2) - DX(2)
        If X(1) = X1(1) And X(2) = X1(2) Then
            Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
            Exit Sub
        End If
        X(1) = X1(1)
        X(2) = X1(2)
        CDX = CABS(DX(1), DX(2))
        DXOLD = CDX
        If Not POLISH% Then
            If CDX <= EPS * CABS(X(1), X(2)) Then
                Erase X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO
                Exit Sub
            End If
        End If
    Next ITER
    Print "too many iterations"
End Sub
Function CABS(A1, A2)
    X = Abs(A1)
    Y = Abs(A2)
    If X = 0 Then
        CABS = Y
    ElseIf Y = 0 Then
        CABS = X
    ElseIf X > Y Then
        CABS = X * Sqr(1 + Sqr(Y / X))
    Else
        CABS = Y * Sqr(1 + Sqr(X / Y))
    End If
End Function
Function CDIV1(A1, A2, B1, B2)
    If Abs(B1) >= Abs(B2) Then
        R = B2 / B1
        DEN = B1 + R * B2
        CDIV1 = (A1 + A2 * R) / DEN
    Else
        R = B1 / B2
        DEN = B2 + R * B1
        CDIV1 = (A1 * R + A2) / DEN
    End If
End Function
Function CDIV2(A1, A2, B1, B2)
    If Abs(B1) >= Abs(B2) Then
        R = B2 / B1
        DEN = B1 + R * B2
        CDIV2 = (A2 - A1 * R) / DEN
    Else
        R = B1 / B2
        DEN = B2 + R * B1
        CDIV2 = (A2 * R - A1) / DEN
    End If
End Function
Function CSQR1(X, Y)
    If X = 0 And Y = 0 Then
        U = 0
    Else
        If Abs(X) >= Abs(Y) Then
            W = Sqr(Abs(X)) * Sqr(0.5 * (1 + Sqr(1 + Sqr(Abs(Y / X)))))
        Else
            R = Abs(X / Y)
            W = Sqr(Abs(Y)) * Sqr(0.5 * (R + Sqr(1 + Sqr(R))))
        End If
        If X >= 0 Then
            U = W
            V = Y / (2 * U)
        Else
            If Y >= 0 Then
            V = W
        Else
            V = -W
        End If
        U = Y / (2 * V)
        End If
    End If
    CSQR1 = U
End Function
Function CSQR2(X, Y)
    If X = 0 And Y = 0 Then
        V = 0
    Else
        If Abs(X) >= Abs(Y) Then
            W = Sqr(Abs(X)) * Sqr(0.5 * (1 + Sqr(1 + Sqr(Abs(Y / X)))))
        Else
            R = Abs(X / Y)
            W = Sqr(Abs(Y)) * Sqr(0.5 * (R + Sqr(1 + Sqr(R))))
        End If
        If X >= 0 Then
            U = W
            V = Y / (2 * U)
        Else
            If Y >= 0 Then
            V = W
        Else
            V = -W
        End If
        U = Y / (2 * V)
        End If
    End If
    CSQR2 = V
End Function