interpmodule.bas

几种常用的内插数值计算方法

'          h     - Integer型变量，等距结点的步长
'          x0    - Double型变量，存放等距n个结点中第一个结点的值
'          y     - Double型一维数组，长度为n，存放给定的n个等距结点的函数值y(i)，y(i) = f(x(i)), i=1,2,...,n
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INEdPq(n As Integer, h As Double, x0 As Double, y() As Double, t As Double) As Double
    Dim i As Integer, j As Integer, k As Integer, m As Integer, l As Integer
    Dim z As Double, hh As Double, xi As Double, xj As Double, b(8) As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INEdPq = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INEdPq = z
        Exit Function
    End If
        
    ' 开始插值
    If (n <= 8) Then
        k = 0
        m = n
    Else
        If (t < (x0 + 4# * h)) Then
            k = 0
            m = 8
        Else
            If (t > (x0 + (n - 5) * h)) Then
                k = n - 8
                m = 8
            Else
                k = Int((t - x0) / h) - 3
                m = 8
            End If
        End If
    End If
        
    b(1) = y(k + 1)
    For i = 2 To m
        hh = y(i + k)
        l = 0
        j = 1
        
        While ((l = 0) And (j <= i - 1))
            If (Abs(hh - b(j)) + 1# = 1#) Then
                l = 1
            Else
                xi = x0 + (i + k - 1) * h
                xj = x0 + (j + k - 1) * h
                hh = (xi - xj) / (hh - b(j))
            End If
            
            j = j + 1
        Wend
        
        b(i) = hh
        If (l <> 0) Then b(i) = 1E+35
    Next i
    
    z = b(m)
    For i = m - 1 To 1 Step -1
      z = b(i) + (t - (x0 + (i + k - 1) * h)) / z
    Next i
    
    ' 返回结果
    INEdPq = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INHermite
'  功能：  进行埃尔米特不等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'          y    - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'         dy     - Double型一维数组，长度为n，存放给定的n个结点的一阶导数值y'(i)，y'(i) = f'(x(i)), i=1,2,...,n
'          t     - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INHermite(n As Integer, x() As Double, y() As Double, dy() As Double, t As Double) As Double
    Dim i As Integer, j As Integer
    Dim z As Double, p As Double, q As Double, s As Double
    
    ' 初值
    z = 0#
    
    ' 循环计算
    For i = 1 To n
        s = 1#
        
        For j = 1 To n
          If (j <> i) Then s = s * (t - x(j)) / (x(i) - x(j))
        Next j
        
        s = s * s
        p = 0#
        For j = 1 To n
          If (j <> i) Then p = p + 1# / (x(i) - x(j))
        Next j
        
        q = y(i) + (t - x(i)) * (dy(i) - 2# * y(i) * p)
        z = z + q * s
    Next i
    
   ' 返回结果
   INHermite = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdHermite
'  功能：  进行埃尔米特等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          h     - Integer型变量，等距结点的步长
'          x0    - Double型变量，存放等距n个结点中第一个结点的值
'          y     - Double型一维数组，长度为n，存放给定的n个等距结点的函数值y(i)，y(i) = f(x(i)), i=1,2,...,n
'         dy     - Double型一维数组，长度为n，存放给定的n个结点的一阶导数值y'(i)，y'(i) = f'(x(i)), i=1,2,...,n
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INEdHermite(n As Integer, h As Double, x0 As Double, y() As Double, dy() As Double, t As Double) As Double
    Dim i As Integer, j As Integer
    Dim z As Double, s As Double, p As Double, q As Double
    
    ' 初值
    z = 0#
    
    ' 循环计算
    For i = 1 To n
        s = 1#
        q = x0 + (i - 1) * h
        
        For j = 1 To n
            p = x0 + (j - 1) * h
            If (j <> i) Then s = s * (t - p) / (q - p)
        Next j
        
        s = s * s
        p = 0#
        
        For j = 1 To n
          If (j <> i) Then p = p + 1# / (q - (x0 + (j - 1) * h))
        Next j
        
        q = y(i) + (t - q) * (dy(i) - 2# * y(i) * p)
        z = z + q * s
    Next i
    
   ' 返回结果
    INEdHermite = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INAitken
'  功能：  进行埃特金不等距逐步插值
'  参数：  n     - Integer型变量，给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'          y    - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'          t     - Double型变量，存放指定的插值点的值
'         eps     - Double型变量，精度控制参数
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INAitken(n As Integer, x() As Double, y() As Double, t As Double, eps As Double) As Double
    Dim i As Integer, j As Integer, k As Integer, m As Integer, l As Integer
    Dim z As Double, xx(10) As Double, yy(10) As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INAitken = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INAitken = z
        Exit Function
    End If
    
    ' 开始插值
    m = 10
    If (m > n) Then m = n
    
    If (t <= x(1)) Then
        k = 1
    Else
        If (t >= x(n)) Then
            k = n
        Else
            k = 1
            j = n
        
            While ((k - j <> 1) And (k - j <> -1))
                l = (k + j) / 2
                If (t < x(l)) Then
                    j = l
                Else
                    k = l
                End If
            Wend
            
            If (Abs(t - x(l)) > Abs(t - x(j))) Then k = j
        End If
    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
        
        xx(i) = x(k)
        yy(i) = y(k)
        l = l + 1
        j = -j
    Next i
    
    i = 0
    Do
        i = i + 1
        z = yy(i + 1)
        For j = 1 To i
          z = yy(j) + (t - xx(j)) * (yy(j) - z) / (xx(j) - xx(i + 1))
        Next j
        
        yy(i + 1) = z
    Loop While ((i <> m - 1) And (Abs(yy(i + 1) - yy(i)) > eps))
    
    ' 返回结果
    INAitken = z

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdAitken
'  功能：  进行埃特金等距逐步插值
'  参数：  n     - Integer型变量，给定结点的点数
'          h     - Integer型变量，等距结点的步长
'          x0    - Double型变量，存放等距n个结点中第一个结点的值
'          y     - Double型一维数组，长度为n，存放给定的n个等距结点的函数值y(i)，y(i) = f(x(i)), i=1,2,...,n
'         dy     - Double型一维数组，长度为n，存放给定的n个结点的一阶导数值y'(i)，y'(i) = f'(x(i)), i=1,2,...,n
'          t   - Double型变量，存放指定的插值点的值
'         eps     - Double型变量，精度控制参数
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INEdAitken(n As Integer, h As Double, x0 As Double, y() As Double, t As Double, eps As Double) As Double
    Dim i As Integer, j As Integer, k As Integer, m As Integer, l As Integer
    Dim z As Double, xx(10) As Double, yy(10) As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INEdAitken = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INEdAitken = z
        Exit Function
    End If
    
    ' 开始插值
    m = 10
    If (m > n) Then m = n
    
    If (t <= x0) Then
        k = 1
    Else
        If (t >= x0 + (n - 1) * h) Then
            k = n
        Else
            k = 1
            j = n
        
            While ((k - j <> 1) And (k - j <> -1))
                l = (k + j) / 2
                If (t < x0 + (l - 1) * h) Then
                    j = l
                Else
                    k = l
                End If
            Wend
            
            If (Abs(t - (x0 + (l - 1) * h)) > Abs(t - (x0 + (j - 1) * h))) Then k = j
            
        End If
    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
        
        xx(i) = x0 + (k - 1) * h
        yy(i) = y(k)
        l = l + 1
        j = -j
    Next i
    
    i = 0
    
    Do
        i = i + 1
        z = yy(i + 1)
        
        For j = 1 To i
            z = yy(j) + (t - xx(j)) * (yy(j) - z) / (xx(j) - xx(i + 1))
        Next j
        
        yy(i + 1) = z
    Loop While ((i <> m - 1) And (Abs(yy(i + 1) - yy(i)) > eps))
    
    ' 返回结果
    INEdAitken = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INAkima
'  功能：  光滑不等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'          y     - Double型一维数组，长度为n，存放给定的n个结点的函数值y(i)，y(i) = f(x(i)), i=1,2,...,n
'          k     - Integer型变量，控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数
'                 s1，s2，s3，s4；若k<0，则需要计算指定插值点t处的函数近似值f(t)，并计算所在子区间的三
'                 次多项式系数s1，s2，s3，s4
'          t     - Double型变量，存放指定的插值点的值,若k>=0，此参数不起作用，可为任意值
'          s     - Double型一维数组，长度为5，其中s1，s2，s3，s4返回三次多项式的系数，s5返回指定插值点t处的
'                 函数近似值f(t)（k<0时）或任意值（k>=0时）
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub INAkima(n As Integer, x() As Double, y() As Double, k As Integer, t As Double, s() As Double)
    Dim kk As Integer, l As Integer, m As Integer
    Dim u(5) As Double, p As Double, q As Double
    
    ' 初值
    s(5) = 0#
    s(1) = 0#
    s(2) = 0#
    s(3) = 0#
    s(4) = 0#
    
    ' 特例处理
    If (n < 1) Then
        Exit Sub
    End If
    
    If (n = 1) Then
        s(1) = y(1)
        s(5) = y(1)
        Exit Sub
    End If
    
    If (n = 2) Then
        s(1) = y(1)
        s(2) = (y(2) - y(1)) / (x(2) - x(1))
        If (k < 0) Then
          s(5) = (y(1) * (t - x(2)) - y(2) * (t - x(1))) / (x(1) - x(2))
        End If
        Exit Sub
    End If
    
    ' 开始插值
    If (k < 0) Then
        If (t <= x(2)) Then
            kk = 0
        Else
            If (t >= x(n)) Then