d10r11.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5595
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   5595
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3000
      TabIndex        =   0
      Top             =   4920
      Width           =   1335
   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 D10R11
    'Driver for routine ZROOTS
    M = 4
    M1 = M + 1
    Dim A(2, 5), X, ROOTS(2, 4)
    For J = 1 To M1
        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
    POLISH% = 0
    Call ZROOTS(A(), M, ROOTS(), POLISH%)
    Print Tab(5); "Unpolished roots:"
    Print Tab(3); "Root #        Real          Imag."
    For I = 1 To M
        Print Tab(5); Format$(I, "##");
        Print Tab(15); Format$(ROOTS(1, I), "#.####00");
        Print Tab(29); Format$(ROOTS(2, I), "#.####00")
    Next I
    Print Tab(5)
    Print Tab(5); "Corrupted roots:"
    For I = 1 To M
        ROOTS(1, I) = ROOTS(1, I) * (1# + 0.01 * I)
        ROOTS(2, I) = ROOTS(2, I) * (1# + 0.01 * I)
    Next
    Print Tab(3); "Roots #       Real          Imag."
    For I = 1 To M
        Print Tab(5); Format$(I, "##");
        Print Tab(15); Format$(ROOTS(1, I), "#.#####0");
        Print Tab(29); Format$(ROOTS(2, I), "#.#####0")
    Next I
    POLISH% = -1
    Call ZROOTS(A(), M, ROOTS(), POLISH%)
    Print Tab(5)
    Print Tab(5); "Polished roots:"
    Print Tab(3); "Roots #       Real          Imag."
    For I = 1 To M
        Print Tab(5); Format$(I, "##");
        Print Tab(15); Format$(ROOTS(1, I), "#.#####0");
        Print Tab(29); Format$(ROOTS(2, I), "#.#####0")
    Next I
End Sub
Sub ZROOTS(A(), M, ROOTS(), POLISH%)
    EPS = 0.000001
    Dim AD(2, 101), X(2), B(2), C(2), DUM(2)
    For J = 1 To M + 1
        AD(1, J) = A(1, J)
        AD(2, J) = A(2, J)
    Next J
    For J = M To 1 Step -1
        X(1) = 0
        X(2) = 0
        Call LAGUER(AD(), J, X(), EPS, 0)
        If Abs(X(2)) <= 2# * EPS ^ 2 * Abs(X(1)) Then X(2) = 0
        ROOTS(1, J) = X(1)
        ROOTS(2, J) = X(2)
        B(1) = AD(1, J + 1)
        B(2) = AD(2, J + 1)
        For JJ = J To 1 Step -1
            C(1) = AD(1, JJ)
            C(2) = AD(2, JJ)
            AD(1, JJ) = B(1)
            AD(2, JJ) = B(2)
            DUM1 = B(1)
            B(1) = X(1) * DUM1 - X(2) * B(2) + C(1)
            B(2) = X(2) * DUM1 + X(1) * B(2) + C(2)
        Next JJ
    Next J
    If POLISH% Then
        For J = 1 To M
            DUM(1) = ROOTS(1, J)
            DUM(2) = ROOTS(2, J)
            Call LAGUER(A(), M, DUM(), EPS, -1)
        Next J
    End If
    For J = 2 To M
        X(1) = ROOTS(1, J)
        X(2) = ROOTS(2, J)
        For I = J - 1 To 1 Step -1
            If ROOTS(1, I) <= X(1) Then Exit For
            ROOTS(1, I + 1) = ROOTS(1, I)
            ROOTS(2, I + 1) = ROOTS(2, I)
        Next I
        If ROOTS(1, I) > X(1) Then I = 0
        ROOTS(1, I + 1) = X(1)
        ROOTS(2, I + 1) = X(2)
    Next J
    Erase DUM, C, B, X, AD
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