曲线_插值m2.bas

<VB数理统计实用算法>书中的算法源程序

Attribute VB_Name = "modMethod"
'曲线_插值
Option Explicit

'线性插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub LIP(X() As Double, Y() As Double, T As Double, F As Double)
    Dim I As Integer, N As Integer
    On Error GoTo errL
    N = UBound(X, 1)                        '数据点数
    For I = 1 To N - 1
        If T < X(I + 1) Then
            F = Y(I) + (Y(I + 1) - Y(I)) / (X(I + 1) - X(I)) * (T - X(I))
            Exit Sub
        End If
        If I = N - 1 Then
            F = Y(N - 1) + (Y(N) - Y(N - 1)) / (X(N) - X(N - 1)) * (T - X(N - 1))
        End If
    Next I
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'一元三点插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub QIP(X() As Double, Y() As Double, T As Double, F As Double)
    Dim I As Integer, N As Integer
    Dim U As Double, V As Double, W As Double
    Dim X0 As Double, X1 As Double, X2 As Double
    On Error GoTo errL
    N = UBound(X, 1)                        '数据点数
    For I = 1 To N - 2
        If T < X(I + 1) Then
            X0 = X(I): X1 = X(I + 1): X2 = X(I + 2)
            U = (T - X1) * (T - X2) / ((X0 - X1) * (X0 - X2))
            V = (T - X0) * (T - X2) / ((X1 - X0) * (X1 - X2))
            W = (T - X0) * (T - X1) / ((X2 - X0) * (X2 - X1))
            F = U * Y(I) + V * Y(I + 1) + W * Y(I + 2)
            Exit Sub
        End If
        If I = N - 2 Then
            X0 = X(N - 2): X1 = X(N - 1): X2 = X(N)
            U = (T - X1) * (T - X2) / ((X0 - X1) * (X0 - X2))
            V = (T - X0) * (T - X2) / ((X1 - X0) * (X1 - X2))
            W = (T - X0) * (T - X1) / ((X2 - X0) * (X2 - X1))
            F = U * Y(N - 2) + V * Y(N - 1) + W * Y(N)
        End If
    Next I
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'拉格朗日插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub LAGRAN(X() As Double, Y() As Double, T As Double, F As Double)
    Dim I As Integer, J As Integer, N As Integer
    Dim P As Double
    On Error GoTo errL
    N = UBound(X, 1)                        '数据点数
    F = 0#
    For I = 1 To N
        P = 1#
        For J = 1 To N
            If (J <> I) Then P = P * (T - X(J)) / (X(I) - X(J))
        Next J
        F = F + P * Y(I)
    Next I
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'牛顿插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub NEWTON(X() As Double, Y() As Double, T As Double, F As Double)
    Dim I As Integer, J As Integer, N As Integer
    Dim B(1000) As Double, D(1000) As Double
    On Error GoTo errL
    N = UBound(X, 1)                        '数据点数
    For I = 1 To N
        D(I) = Y(I)
    Next I
    For I = 1 To N - 1
        For J = N To I + 1 Step -1
            D(J) = (D(J - 1) - D(J)) / (X(J - I) - X(J))
        Next J
    Next I
    B(N) = D(N)
    For I = N - 1 To 1 Step -1
        B(I) = D(I) + (T - X(I)) * B(I + 1)
    Next I
    F = B(1)
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'埃特金插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub AITKEN(X() As Double, Y() As Double, T As Double, F As Double)
    Dim XM(1 To 1000) As Double, YM(1 To 1000) As Double
    Dim I As Integer, J As Integer, K As Integer
    Dim L As Integer, N As Integer, M As Integer
    On Error GoTo errL
    N = UBound(X, 1)        '数据点数
    M = 10                  '用最靠近插值点T的M个数据点作埃特金插值
    If M > N Then M = N
    If T <= X(1) Then
        K = 1
    ElseIf T > X(N) Then
        K = N
    Else
        K = 1: J = N
10:
        If Abs(K - J) <> 1 Then
            L = (K + J) / 2
            If T < X(L) Then J = L Else K = L
            GoTo 10
        End If
        If Abs(T - X(L)) > Abs(T - X(J)) Then K = J
    End If
    J = 1: L = 0
    For I = 1 To M
        K = K + J * L
        If K < 1 Or K > N Then
            L = L + 1: J = -J: K = K + J * L
        End If
        XM(I) = X(K): YM(I) = Y(K)
        L = L + 1: J = -J
    Next I
    For I = 2 To M
        F = YM(I)
        For J = 2 To I
            F = YM(J - 1) + (T - XM(J - 1)) * _
                (YM(J - 1) - F) / (XM(J - 1) - XM(I))
        Next J
        YM(I) = F
    Next I
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'三次样条函数插值
'X：数据点数组
'Y：函数值数组
'T：插值点
'F：插值点函数值
Public Sub SPLINE(X() As Double, Y() As Double, T As Double, F As Double)
    Dim I As Integer, J As Integer, K As Double
    Dim N As Integer, N1 As Integer, N2 As Integer
    Dim Z As Double, H1 As Double, H2 As Double, H3 As Double, H4 As Double
    Dim H(1000) As Double, DY(1000) As Double, F2(1000) As Double
    Dim S(1000) As Double, E(1000) As Double
    On Error GoTo errL
    N = UBound(X, 1)                        '数据点数
    N1 = N - 1: N2 = N - 2
    F2(1) = 0: F2(N) = 0
    For I = 1 To N1
        H(I) = X(I + 1) - X(I)
        DY(I) = (Y(I + 1) - Y(I)) / H(I)
    Next I
    For I = 2 To N1
        F2(I) = 6# * (DY(I) - DY(I - 1))
    Next I
    Z = 0.5 / (H(1) + H(2)): S(1) = -H(2) * Z
    E(1) = F2(2) * Z: K = 1
    For I = 2 To N2
        J = I + 1
        Z = 1# / (2# * (H(I) + H(J)) + H(I) * S(K))
        S(I) = -H(J) * Z
        E(I) = (F2(J) - H(I) * E(K)) * Z
        K = 1
    Next I
    F2(N1) = E(N2)
    For I = N2 To 2 Step -1
        K = I - 1: F2(I) = S(K) * F2(I + 1) + E(K)
    Next I
    For I = 1 To N1
        S(I) = (F2(I + 1) - F2(2)) / H(I)
    Next I
    I = 2: K = 1
LL:
    If I < N And T > X(I) Then K = I: I = I + 1: GoTo LL
    If I = N Or T <= X(I) Then
        H1 = T - X(K): H2 = T - X(I)
        H3 = H1 * H2: H4 = F2(K) + H1 * S(K)
        Z = (F2(I) + F2(K) + H4) / 6#
        F = Y(K) + H1 * DY(K) + H3 * Z
    End If
    Exit Sub
errL:
    MsgBox "不同的数据点有相同的X坐标，造成除数为0"
End Sub

'计算等距插值点的函数值（等距化时使用）
'X：数组，观测数据的X坐标
'Y：数组，观测数据的Y坐标
'B：数组，保存等距多点的函数值
'LLL：方法
Public Sub Equal(X() As Double, Y() As Double, B() As Double, LLL As Integer)
    Dim N As Integer, M As Integer, I As Integer
    Dim T As Double, F As Double
    Dim miX As Double, maX As Double, DX As Double
    N = UBound(X, 1)                        'N为观测点个数
    M = UBound(B, 2)                        '网格的行数
    miX = X(1)                              'X坐标最小值
    maX = X(N)                              'X坐标最大值
    DX = (maX - miX) / (M - 1)              '网格在X方向上的增量
    For I = 1 To M
        T = miX + DX * (I - 1)
        Select Case LLL
            Case 0
                LIP X, Y, T, F              '线性插值
            Case 1
                QIP X, Y, T, F              '一元三点插值
            Case 2
                LAGRAN X, Y, T, F           '拉格朗日插值
            Case 3
                NEWTON X, Y, T, F           '牛顿插值
            Case 4
                AITKEN X, Y, T, F           '埃特金插值
            Case 5
                SPLINE X, Y, T, F           '三次样条函数插值
        End Select
        B(1, I) = T                         '插值点X坐标
        B(2, I) = F                         '插值点函数值
    Next I
End Sub

'单点插值
'X：数组，观测数据的X坐标
'Y：数组，观测数据的Y坐标
'T：插值点X坐标
'F：单点插值的函数值
'LLL：方法
Public Sub OneP(X() As Double, Y() As Double, _
        T As Double, F As Double, LLL As Integer)
    Select Case LLL
        Case 0
            LIP X, Y, T, F                  '线性插值
        Case 1
            QIP X, Y, T, F                  '一元三点插值
        Case 2
            LAGRAN X, Y, T, F               '拉格朗日插值
        Case 3
            NEWTON X, Y, T, F               '牛顿插值
        Case 4
            AITKEN X, Y, T, F               '埃特金插值
        Case 5
            SPLINE X, Y, T, F               '三次样条函数插值
    End Select
End Sub