libmathfuncs.bas

本系统是给大庆油田做的一个示例程序

    Dim X1(5) As Single, Y1(5) As Single
    Dim A As Single, B As Single
    Dim yy(5000) As Single
    
    
    yy(1) = Y(1)
    yy(2) = Y(2)
    yy(M - 1) = Y(M - 1)
    yy(M) = Y(M)
    
    For j = 3 To M - 2
        k = 5
        For i = 1 To 5
            X1(i) = x(j - 3 + i)
            Y1(i) = Y(j - 3 + i)
        Next i
        Call Regression_Linear(k, X1, Y1, A, B)
        yy(j) = A + B * x(j)
    Next j
    For j = 1 To M
        Y(j) = yy(j)
    Next j
End Sub


Rem  线性回归系数
   '文件名： 一元线性元回归  Y=A+Bx
   '函数名： Regression_Linear( K , X(), Y(), A())
   '参  数：K —— 整型变量，输入参数。观测点数。
   '        X() —— 实型1维数组，输入参数。自变量因子的k次观测值
   '        Y() —— 实型1维数组，输入参数。因变量因子的k次观测值
Sub Regression_Linear(k As Integer, x() As Single, Y() As Single, A As Single, B As Single)
    Dim j As Integer
    Dim Xm As Single, Ym As Single
    Dim Sum1 As Single, Sum2 As Single
    
    Sum1 = 0
    Sum2 = 0
    For j = 1 To k
        Sum1 = Sum1 + x(j)
        Sum2 = Sum2 + Y(j)
    Next j
    
    Xm = Sum1 / k
    Ym = Sum2 / k
    
    Sum1 = 0
    Sum2 = 0
    For j = 1 To k
        Sum1 = Sum1 + (x(j) - Xm) * (Y(j) - Ym)
        Sum2 = Sum2 + (x(j) - Xm) ^ 2
    Next j
    
    B = Sum1 / Sum2
    A = Ym - B * Xm
End Sub
Rem 斜率小于给定值区间载荷位移曲线的抛物线光滑处理
Sub PumpCard_ParabolaSmoothing(Ncal As Integer, Kmin As Single, UnitXPump() As Single, UnitPPump() As Single)
    
    Dim j As Integer, k As Integer, i As Integer
    Dim iPoint_KLessKmin(600) As Integer
    Dim Ncal_KLessKmin As Integer, NRegion As Integer
    Dim iPoint_Start(600) As Integer, iPoint_End(600) As Integer
    
    Dim A(5) As Single
    Dim M As Integer, n As Integer
    Dim x(1500) As Single, Y(1500) As Single, T(1500) As Single
    Dim F(1500) As Single, F1(1500) As Single, F2(1500) As Single
    
    n = Ncal
    M = n
    For j = 1 To n
        x(j) = UnitXPump(j)
        Y(j) = UnitPPump(j)
    Next j
    
    For j = 1 To M
         T(j) = 1 / (M - 1) * (j - 1)
    Next j
    
    Call ChZh(n, x, Y, M, T, F, F1, F2)
   
    For j = 1 To M
        UnitXPump(j) = T(j)
        UnitPPump(j) = F(j)
        If F1(j) >= 4.5 Then F1(j) = 4.5
        If F1(j) <= -4.5 Then F1(j) = -4.5
        If F2(j) >= 4.5 Then F2(j) = 4.5
        If F2(j) <= -4.5 Then F2(j) = -4.5
    Next j
    
    k = 0
    For j = 1 To Ncal
        If Abs(F1(j)) <= Kmin Then
            k = k + 1
            iPoint_KLessKmin(k) = j
        End If
    Next j
    Ncal_KLessKmin = k
    
    iPoint_Start(1) = iPoint_KLessKmin(1)
    k = 1
    For j = 2 To Ncal_KLessKmin
        If (iPoint_KLessKmin(j) - iPoint_KLessKmin(j - 1)) > 1 Then
            k = k + 1
            iPoint_Start(k) = iPoint_KLessKmin(j)
            iPoint_End(k - 1) = iPoint_KLessKmin(j - 1)
        Else
        End If
    Next j
    iPoint_End(k) = iPoint_KLessKmin(Ncal_KLessKmin)
    NRegion = k
    
    For j = 1 To NRegion
        
        If iPoint_Start(j) = iPoint_End(j) Then
            UnitPPump(iPoint_Start(j)) = (UnitPPump(iPoint_Start(j) - 1) + UnitPPump(iPoint_Start(j) + 1)) / 2
        End If
        
        If iPoint_End(j) - iPoint_Start(j) = 1 Then
            If iPoint_Start(j) = 1 Then
                UnitPPump(2) = (UnitPPump(1) + UnitPPump(3)) / 2
            End If
            If iPoint_Start(j) = Ncal Then
                UnitPPump(Ncal - 1) = (UnitPPump(Ncal) + UnitPPump(Ncal - 2)) / 2
            End If
            If iPoint_Start(j) > 1 And iPoint_End(j) < Ncal Then
                UnitPPump(iPoint_Start(j)) = UnitPPump(iPoint_Start(j) - 1) + _
                           (UnitPPump(iPoint_End(j) + 1) - UnitPPump(iPoint_Start(j) - 1)) / 3 * 1
                UnitPPump(iPoint_End(j)) = UnitPPump(iPoint_Start(j) - 1) + _
                           (UnitPPump(iPoint_End(j) + 1) - UnitPPump(iPoint_Start(j) - 1)) / 3 * 1
            End If
            
        End If
        
        If iPoint_End(j) - iPoint_Start(j) >= 2 Then
            k = iPoint_End(j) - iPoint_Start(j) + 1
            For i = 1 To k
                x(i) = UnitXPump(iPoint_Start(j) + i - 1)
                F(i) = UnitPPump(iPoint_Start(j) + i - 1)
            Next i
            Call Regression_Parabola(k, x, F, A)
            For i = 1 To k
                UnitPPump(iPoint_Start(j) + i - 1) = A(0) + A(1) * UnitXPump(iPoint_Start(j) + i - 1) + _
                       A(2) * UnitXPump(iPoint_Start(j) + i - 1) ^ 2
            Next i
        End If
    Next j
End Sub

Sub Gauss(d(), e(), F())
    Dim i As Integer, j As Integer, L As Integer, k As Integer
    Dim N1 As Integer
    
    Dim D11 As Single, aaa1 As Single, s As Single
     
    N1 = UBound(d, 1)
    ReDim A(N1, N1 + 1)
    
    For i = 1 To N1
        A(i, N1 + 1) = e(i, 1)
        For j = 1 To N1
            A(i, j) = d(i, j)
        Next j
    Next i
    For k = 1 To N1 - 1
        D11 = A(k, k)
        L = k
        
        For i = k + 1 To N1
            If Abs(A(i, k)) > Abs(D11) Then
                D11 = A(i, k)
                L = i
            Else
            End If
        Next i
        If D11 = 0 Then
            Stop
        Else
        End If
        If L <> k Then
            For j = k To N1 + 1
                aaa1 = A(L, j)
                A(L, j) = A(k, j)
                A(k, j) = aaa1
            Next j
        Else
        End If
        For j = k + 1 To N1 + 1
            A(k, j) = A(k, j) / A(k, k)
        Next j
        For j = k + 1 To N1 + 1
            For i = k + 1 To N1
                A(i, j) = A(i, j) - A(i, k) * A(k, j)
            Next i
        Next j
    Next k
    F(N1, 1) = A(N1, N1 + 1) / A(N1, N1)
    For i = N1 - 1 To 1 Step -1
        s = 0
        For j = i + 1 To N1
            s = s + A(i, j) * F(j, 1)
        Next j
        F(i, 1) = A(i, N1 + 1) - s
    Next i
End Sub
Sub YtChZh(mpn0, cal0(), cpn, scal0())
    
    Dim i As Integer, j As Integer
    
    Dim tt As Single, hcal0 As Single, xhc As Single, LL  As Single
    Dim xl1 As Single, xl2 As Single
    
    ReDim lamda(0 To mpn0), mu_l(0 To mpn0), dcal0(0 To mpn0)
    ReDim mcal0(0 To mpn0), pcal0(0 To mpn0), q(0 To mpn0)
    
    tt = 100
    hcal0 = tt / mpn0
    For i = 0 To mpn0
        lamda(i) = 0.5
        mu_l(i) = 0.5
    Next i
  
    cal0(0) = cal0(mpn0)
    cal0(mpn0 + 1) = cal0(1)
   
    For i = 1 To mpn0
        dcal0(i) = 3 * (cal0(i + 1) + cal0(i - 1) - 2 * cal0(i)) / hcal0 ^ 2
    Next i
    For i = 1 To mpn0
        pcal0(i) = (dcal0(i) - mu_l(i) * pcal0(i - 1)) / (2 - mu_l(i) * q(i - 1))
        q(i) = lamda(i) / (2 - mu_l(i) * q(i - 1))
    Next i
    mcal0(mpn0) = pcal0(mpn0)
    
    For i = mpn0 - 1 To 0 Step -1
        mcal0(i) = pcal0(i) - q(i) * mcal0(i + 1)
    Next i
    
    xhc = tt / cpn
    For i = 0 To cpn
        LL = xhc * i
        For j = 1 To mpn0
            If LL <= j * hcal0 And LL >= (j - 1) * hcal0 Then
               xl1 = j * hcal0 - LL
               xl2 = LL - (j - 1) * hcal0
               scal0(i) = (mcal0(j - 1) * xl1 ^ 3 + mcal0(j) * xl2 ^ 3) / 6 / hcal0
               scal0(i) = scal0(i) + (cal0(j - 1) - mcal0(j - 1) * hcal0 ^ 2 / 6) * xl1 / hcal0
               scal0(i) = scal0(i) + (cal0(j) - mcal0(j) * hcal0 ^ 2 / 6) * xl2 / hcal0
            Else
            End If
        Next j
    Next i
End Sub
Rem 离散点付氏级数变换
Sub Xn_Fain_DFT(Ncal As Integer, UnitTime() As Single, UnitPPump() As Single, _
                 NcalXn As Integer, An() As Single, Bn() As Single, Xn() As Single, Fain() As Single)
    
    Dim j As Integer, k As Integer
    Dim sum As Single
    
    For k = 0 To NcalXn
        sum = 0
        For j = 1 To Ncal
            sum = sum + UnitPPump(j) * Cos(2 * k * pi * j / Ncal)
        Next j
        An(k) = 2 / Ncal * sum
    Next k
    
    For k = 1 To NcalXn
        sum = 0
        For j = 1 To Ncal
            sum = sum + UnitPPump(j) * Sin(2 * k * pi * j / Ncal)
        Next j
        Bn(k) = 2 / Ncal * sum
    Next k
    
    Xn(0) = An(0)
    For k = 1 To NcalXn
        Xn(k) = Sqr(An(k) ^ 2 + Bn(k) ^ 2)
        Fain(k) = Atn(Bn(k) / An(k))
    Next k
    
    For k = 0 To NcalXn
        Xn(k) = Int(Xn(k) * 1000) / 1000
        Fain(k) = Int(Fain(k) * 180 / pi * 100) / 100
    Next k
End Sub