module1.vb

基于vb.net的最小二乘法多项式系数计算类. 很好用,比Excel还精确.

            p = 0
            For j = 2 To Degree + 1
                p = p + b(j) * (bxOrig(j) - covarOrig(j, 1) * bxOrig(1) / ndata)
            Next j
            q = Sumy2 - bxOrig(1) ^ 2 / ndata
            z = q - p
            v = ndata - Degree - 1
            r2 = p / q
            If v > 0 Then StdErr = (z / v) ^ 0.5 Else StdErr = 0
            covsrt(covar, NumTerms, FitTerm, mfit)
        End Sub

        Public Sub funcs(ByVal x As Double, ByRef afunc() As Double, ByVal NumTerms As Integer)
            'generate array afunc(1) = 1, afunc(2) = x, afunc(3) = x^2, afunc(4)= x^3,...
            ' evalute on  x
            Dim i As Integer
            afunc(1) = 1.0
            For i = 2 To NumTerms
                afunc(i) = x * afunc(i - 1)
            Next i
        End Sub

        Sub covsrt(ByVal covar(,) As Double, ByVal NumTerms As Integer, _
                   ByVal FitTerm() As Integer, ByVal mfit As Integer)
            'Expand in storage the covariance matrix covar, so as to take into account parameters that
            'are being held fixed. (For the latter, return zero covariances.)

            Dim i, j, k, h As Integer

            For i = mfit + 1 To NumTerms
                For j = 1 To i
                    covar(i, j) = 0.0
                    covar(j, i) = 0.0
                Next j
            Next i
            k = mfit
            For h = 1 To NumTerms
                j = NumTerms - h + 1
                If FitTerm(j) Then
                    For i = 1 To NumTerms
                        swap(covar(i, k), covar(i, j))
                    Next i
                    For i = 1 To NumTerms
                        swap(covar(k, i), covar(j, i))
                    Next i
                    k = k - 1
                End If
            Next h
        End Sub
        Public Sub swap(ByRef a As Double, ByRef b As Double)
            Dim tmp As Double
            tmp = a
            a = b
            b = tmp
        End Sub
        Public Sub EvaluateX()
            Dim x, y As Double
            Dim xStr As String
            Dim EvalVec(20) As Double
            Dim j As Integer

            xStr = InputBox("Interpolation (enter x value)")
            x = CDbl(xStr)
            funcs(x, EvalVec, Degree + 1)
            y = 0
            For j = 1 To Degree + 1
                y = y + b(j) * EvalVec(j)
            Next j
            MsgBox("For x= " & Format(x, "0.000") & vbCrLf & "y= " & Format(y, "0.000"))
        End Sub
        Public Sub EvaluateY()
            Dim x, y As Double
            Dim xStr, str As String
            Dim i As Integer
            Dim brtw As New Bairstow

            xStr = InputBox("Interpolation (enter y value)")
            y = CDbl(xStr)

            brtw.degree = Degree
            For i = 0 To Degree + 1
                brtw.poly(i) = b(i + 1)  'poly begin from 0, b(i) from 1
            Next
            brtw.poly(0) = b(1) - y

            brtw.brstow()
            For i = 0 To brtw.NumRoots - 1
                If brtw.roots(i, 2) = 0 Then 'only reals roots
                    str = str & "x" & i + 1 & " = " & Format(brtw.roots(i, 1), "0.000") & vbCrLf
                End If

            Next
            MsgBox("For y= " & Format(y, "0.000") & vbCrLf & str)
        End Sub
    End Class

    Public Class GaussJordan
        Private n As Integer
        Private a(50, 50), b(50) As Double

        Public Property aMat() As Array
            'input coefficients of equations
            ' and return inverse matrix when execute gaussj() method 
            Get
                Return a
            End Get
            Set(ByVal value As Array)
                a = value
            End Set
        End Property

        Public Property bVec() As Array
            'input constants of equations
            ' and return solution when execute gaussj() method 
            Get
                Return b
            End Get
            Set(ByVal value As Array)
                b = value
            End Set
        End Property

        Public WriteOnly Property NumEq() As Integer
            ' input number of equations
            Set(ByVal value As Integer)
                n = value
            End Set
        End Property

        Sub gaussj()
            Dim i, icol, irow, j, k, h, f As Integer
            Dim pivinv, big, dum As Double
            Dim indxc(50), indxr(50), ipiv(50) As Integer

            For j = 1 To n
                ipiv(j) = 0
            Next j
            For i = 1 To n
                big = 0.0
                For j = 1 To n
                    If ipiv(j) <> 1 Then
                        For k = 1 To n
                            If ipiv(k) = 0 Then
                                If Abs(a(j, k)) >= big Then
                                    big = Abs(a(j, k))
                                    irow = j
                                    icol = k
                                End If
                            ElseIf ipiv(k) > 1 Then
                                MsgBox("pause 1 in GAUSSJ - singular matrix")
                                Exit Sub
                            End If
                        Next k
                    End If
                Next j
                ipiv(icol) = ipiv(icol) + 1
                If irow <> icol Then
                    For h = 1 To n
                        dum = a(irow, h)
                        a(irow, h) = a(icol, h)
                        a(icol, h) = dum
                    Next h
                    dum = b(irow)
                    b(irow) = b(icol)
                    b(icol) = dum
                End If
                indxr(i) = irow
                indxc(i) = icol
                If a(icol, icol) = 0.0 Then
                    MsgBox("pause 2 in GAUSSJ - singular matrix")
                    Exit Sub
                End If
                pivinv = 1.0 / a(icol, icol)
                a(icol, icol) = 1.0
                For h = 1 To n
                    a(icol, h) = a(icol, h) * pivinv
                Next h
                b(icol) = b(icol) * pivinv
                For f = 1 To n
                    If f <> icol Then
                        dum = a(f, icol)
                        a(f, icol) = 0.0
                        For h = 1 To n
                            a(f, h) = a(f, h) - a(icol, h) * dum
                        Next h
                        b(f) = b(f) - b(icol) * dum
                    End If
                Next f
            Next i
            For i = n To 1 Step -1
                If indxr(i) <> indxc(i) Then
                    For k = 1 To n
                        dum = a(k, indxr(i))
                        a(k, indxr(i)) = a(k, indxc(i))
                        a(k, indxc(i)) = dum
                    Next k
                End If
            Next i
        End Sub
    End Class


    Public Class Bairstow
        Private p(100) As Double '   n+1 polynomial coefficients 
        Private n As Integer     '   degree of polynomial 
        Private tol As Double = 0.000001
        Private nlim As Integer = 50  ' maximum number of iterations 
        Private r, s As Double
        Private rt(20, 2) As Double
        Private count As Integer = 0

        Public Property degree() As Integer
            Get
                Return n
            End Get
            Set(ByVal value As Integer)
                n = value
            End Set
        End Property

        Public Property poly() As Array
            Get
                Return p
            End Get
            Set(ByVal value As Array)
                p = value
            End Set
        End Property
        Public ReadOnly Property roots() As Array
            Get
                Return rt
            End Get

        End Property
        Public ReadOnly Property NumRoots() As Integer
            Get
                Return count
            End Get
        End Property
        Public Property tolerance() As Double
            Get
                Return tol
            End Get
            Set(ByVal value As Double)
                tol = value
            End Set
        End Property
        Public Property iterations() As Integer
            Get
                Return nlim
            End Get
            Set(ByVal value As Integer)
                nlim = value
            End Set
        End Property

        Public Sub brstow()

            Dim i, j, k, np3 As Integer
            Dim a(100), b(100), c(100) As Double

            Dim denom, delr, dels As Double
            Dim discrim, q As Double
            r = 0
            s = 0   '  starting values for quadratic term 

            ' Quadratic factors by Bairstow method
            ' Trial factors: r = 0,  s = 0

            np3 = n + 3
            For i = 0 To n
                a(i) = p(i)
            Next i

            ' begin iteration 
            k = 0
            Do While k <= nlim
                For i = 0 To np3 - 1
                    b(i) = 0.0
                    c(i) = 0.0
                Next i

                ' compute b and c arrays 
                For j = n To 0 Step -1
                    b(j) = a(j) + r * b(j + 1) + s * b(j + 2)
                    c(j) = b(j) + r * c(j + 1) + s * c(j + 2)
                Next j

                ' compute denominator and check if zero. 
                denom = c(2) * c(2) - c(3) * c(1)

                If denom <> 0.0 Then
                    delr = (b(0) * c(3) - b(1) * c(2)) / denom
                    dels = (b(1) * c(1) - b(0) * c(2)) / denom
                    r += delr
                    s += dels
                    If Abs(delr) + Abs(dels) < tol Then
                        'fit in k + 1  iterations
                        'x^2 -  rx -  s
                        check_quadratic(1.0, -r, -s)
                        n -= 2
                        np3 = n + 3
                        For i = 0 To np3 - 1
                            a(i) = b(i + 2)
                        Next i
                        Select Case n
                            Case 0
                                Exit Sub
                            Case 1
                                'a(1)x + a(0)
                                check_quadratic(0.0, a(1), a(0))
                                Exit Sub
                            Case 2
                                ' a(2)x^2 + a(1)x + a(0)
                                check_quadratic(a(2), a(1), a(0))
                                Exit Sub
                            Case Else
                                k = -1 ' reset iteration count
                        End Select
                    End If
                Else
                    ' if denom is zero, perturb r and s values 
                    r += 1.0
                    s += 1.0
                    k = -1
                End If
                k += 1
            Loop
            MsgBox("Bairstow: iteration count of " & nlim & "exceeded")
        End Sub

        '
        '   solve equation a x^2 + b x + c = 0
        '
        Private Sub check_quadratic(ByVal a As Double, ByVal b As Double, ByVal c As Double)
            Dim discrim, q As Double
            Dim re, im As Double
            If Abs(a) < 0.000000000001 Then
                'linear root: -c / b
                rt(count, 1) = -c / b
                count += 1
                Exit Sub
            End If
            ' convert to a x^2 + 2 b x + c = 0 
            b *= 0.5
            discrim = b * b - a * c

            If (discrim > 0.0) Then
                discrim = Sqrt(discrim)
                q = -(b + IIf(b < 0, -discrim, discrim))
                ' real roots: q / a , c / q
                rt(count, 1) = q / a
                count += 1
                rt(count, 1) = c / q
                count += 1
            ElseIf (discrim < 0.0) Then
                re = -b / a
                im = Sqrt(-discrim) / a

                're  + im j
                rt(count, 1) = re
                rt(count, 2) = im
                count += 1
                're - im j
                rt(count, 1) = re
                rt(count, 2) = -im
                count += 1
            Else
                'multiple root: -b / a
                rt(count, 1) = -b / a
                count += 1
            End If
        End Sub

    End Class

End Module