interpmodule.bas

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

                kk = n - 2
            Else
                kk = 1
                m = n
                While (((kk - m) <> 1) And ((kk - m) <> -1))
                    l = (kk + m) / 2
                    If (t < x(l)) Then
                        m = l
                    Else
                        kk = l
                    End If
                Wend
              
                kk = kk - 1
            End If
        End If
    Else
        kk = k
    End If
    
    If (kk >= n - 1) Then kk = n - 2
    
    ' 调用Akima公式
    u(3) = (y(kk + 2) - y(kk + 1)) / (x(kk + 2) - x(kk + 1))
    If (n = 3) Then
        If (kk = 0) Then
            u(4) = (y(3) - y(2)) / (x(3) - x(2))
            u(5) = 2# * u(4) - u(3)
            u(2) = 2# * u(3) - u(4)
            u(1) = 2# * u(2) - u(3)
        Else
            u(2) = (y(2) - y(1)) / (x(2) - x(1))
            u(1) = 2# * u(2) - u(3)
            u(4) = 2# * u(3) - u(2)
            u(5) = 2# * u(4) - u(3)
        End If
    Else
        If (kk <= 1) Then
            u(4) = (y(kk + 3) - y(kk + 2)) / (x(kk + 3) - x(kk + 2))
            If (kk = 1) Then
                u(2) = (y(2) - y(1)) / (x(2) - x(1))
                u(1) = 2# * u(2) - u(3)
                If (n = 4) Then
                    u(5) = 2# * u(4) - u(3)
                Else
                    u(5) = (y(5) - y(4)) / (x(5) - x(4))
                End If
            Else
                u(2) = 2# * u(3) - u(4)
                u(1) = 2# * u(2) - u(3)
                u(5) = (y(4) - y(3)) / (x(4) - x(3))
            End If
        Else
            If (kk >= (n - 3)) Then
                u(2) = (y(kk + 1) - y(kk)) / (x(kk + 1) - x(kk))
                If (kk = (n - 3)) Then
                    u(4) = (y(n) - y(n - 1)) / (x(n) - x(n - 1))
                    u(5) = 2# * u(4) - u(3)
                    If (n = 4) Then
                        u(1) = 2# * u(2) - u(3)
                    Else
                        u(1) = (y(kk) - y(kk - 1)) / (x(kk) - x(kk - 1))
                    End If
                Else
                    u(4) = 2# * u(3) - u(2)
                    u(5) = 2# * u(4) - u(3)
                    u(1) = (y(kk) - y(kk - 1)) / (x(kk) - x(kk - 1))
                End If
            Else
                u(2) = (y(kk + 1) - y(kk)) / (x(kk + 1) - x(kk))
                u(1) = (y(kk) - y(kk - 1)) / (x(kk) - x(kk - 1))
                u(4) = (y(kk + 3) - y(kk + 2)) / (x(kk + 3) - x(kk + 2))
                u(5) = (y(kk + 4) - y(kk + 3)) / (x(kk + 4) - x(kk + 3))
            End If
        End If
    End If
    
    s(1) = Abs(u(4) - u(3))
    s(2) = Abs(u(1) - u(2))
    If ((s(1) + 1# = 1#) And (s(2) + 1# = 1#)) Then
         p = (u(2) + u(3)) / 2#
    Else
        p = (s(1) * u(2) + s(2) * u(3)) / (s(1) + s(2))
    End If
    
    s(1) = Abs(u(4) - u(5))
    s(2) = Abs(u(3) - u(2))
    If ((s(1) + 1# = 1#) And (s(2) + 1# = 1#)) Then
         q = (u(3) + u(4)) / 2#
    Else
        q = (s(1) * u(3) + s(2) * u(4)) / (s(1) + s(2))
    End If
    
    s(1) = y(kk + 1)
    s(2) = p
    s(4) = x(kk + 2) - x(kk + 1)
    s(3) = (3# * u(3) - 2# * p - q) / s(4)
    s(4) = (q + p - 2# * u(3)) / (s(4) * s(4))
    
    If (k < 0) Then
        p = t - x(kk + 1)
        s(5) = s(1) + s(2) * p + s(3) * p * p + s(4) * p * p * p
    End If

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdAkima
'  功能：  光滑等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          h     - Integer型变量，等距结点的步长
'          x0    - Double型变量，存放等距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 INEdAkima(n As Integer, h As Double, x0 As Double, y() As Double, k As Integer, t As Double, s() As Double)
    Dim kk As Integer, m As Integer, l 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(5)
        Exit Sub
    End If
        
    If (n = 2) Then
        s(1) = y(1)
        s(2) = (y(2) - y(1)) / h
        If (k < 0) Then
          s(5) = (y(2) * (t - x0) - y(1) * (t - x0 - h)) / h
        End If
        Exit Sub
    End If
    
    ' 开始插值
    If (k < 0) Then
        If (t <= x0 + h) Then
            kk = 0
        Else
            If (t >= x0 + (n - 1) * h) Then
                kk = n - 2
            Else
                kk = 1
                m = n
                While (((kk - m) <> 1) And ((kk - m) <> -1))
                    l = (kk + m) / 2
                    If (t < x0 + (l - 1) * h) Then
                        m = l
                    Else
                        kk = l
                    End If
                Wend
                
                kk = kk - 1
            End If
        End If
    Else
        kk = k
    End If
    
    If (kk >= n - 1) Then kk = n - 2
    
    ' 调用Akima公式
    u(3) = (y(kk + 2) - y(kk + 1)) / h
    
    If (n = 3) Then
        If (kk = 0) Then
            u(4) = (y(3) - y(2)) / h
            u(5) = 2# * u(4) - u(3)
            u(2) = 2# * u(3) - u(4)
            u(1) = 2# * u(2) - u(3)
        Else
            u(2) = (y(2) - y(1)) / h
            u(1) = 2# * u(2) - u(3)
            u(4) = 2# * u(3) - u(2)
            u(5) = 2# * u(4) - u(3)
        End If
    Else
        If (kk <= 1) Then
            u(4) = (y(kk + 3) - y(kk + 2)) / h
            If (kk = 1) Then
                u(2) = (y(2) - y(1)) / h
                u(1) = 2# * u(2) - u(3)
                If (n = 4) Then
                    u(5) = 2# * u(4) - u(3)
                Else
                    u(5) = (y(5) - y(4)) / h
                End If
            Else
                u(2) = 2# * u(3) - u(4)
                u(1) = 2# * u(2) - u(3)
                u(5) = (y(4) - y(3)) / h
            End If
        Else
            If (kk >= (n - 3)) Then
                u(2) = (y(kk + 1) - y(kk)) / h
                
                If (kk = (n - 3)) Then
                    u(4) = (y(n) - y(n - 1)) / h
                    u(5) = 2# * u(4) - u(3)
                    If (n = 4) Then
                        u(1) = 2# * u(2) - u(3)
                    Else
                        u(1) = (y(kk) - y(kk - 1)) / h
                    End If
                Else
                    u(4) = 2# * u(3) - u(2)
                    u(5) = 2# * u(4) - u(3)
                    u(1) = (y(kk) - y(kk - 1)) / h
                End If
            Else
                u(2) = (y(kk + 1) - y(kk)) / h
                u(1) = (y(kk) - y(kk - 1)) / h
                u(4) = (y(kk + 3) - y(kk + 2)) / h
                u(5) = (y(kk + 4) - y(kk + 3)) / h
            End If
        End If
    End If
    
    s(1) = Abs(u(4) - u(3))
    s(2) = Abs(u(1) - u(2))
    If ((s(1) + 1# = 1#) And (s(2) + 1# = 1#)) Then
         p = (u(2) + u(3)) / 2#
    Else
        p = (s(1) * u(2) + s(2) * u(3)) / (s(1) + s(2))
    End If
    
    s(1) = Abs(u(4) - u(5))
    s(2) = Abs(u(3) - u(2))
    If ((s(1) + 1# = 1#) And (s(2) + 1# = 1#)) Then
         q = (u(3) + u(4)) / 2#
    Else
        q = (s(1) * u(3) + s(2) * u(4)) / (s(1) + s(2))
    End If
    
    s(1) = y(kk + 1)
    s(2) = p
    s(4) = h
    s(3) = (3# * u(3) - 2# * p - q) / s(4)
    s(4) = (q + p - 2# * u(3)) / (s(4) * s(4))
    If (k < 0) Then
        p = t - (x0 + kk * h)
        s(5) = s(1) + s(2) * p + s(3) * p * p + s(4) * p * p * p
    End If
    
End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INSpline1
'  功能：  实现第一种边界条件的三次样条函数插值、微商与积分
'  参数：  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
'          dy    - Double型一维数组，长度为n，调用时，dy(1)存放给定区间的左端点处的一阶导数值，dy(n)存放
'                 给定区间的右端点处的一阶导数值。返回时，存放n个给定点处的一阶导数值y'(i)，i=1,2,…n
'          ddy  - Double型一维数组，长度为n，返回时，存放n个给定点处的二阶导数值y''(i)，i=1,2,…n
'          m    - Integer型变量，指定插值点的个数
'          t     - Double型一维数组，长度为m，存放m个指定的插值点的值。要求x(1)<t(j)<x(n), j=1,2,…,m-1
'          z     - Double型一维数组，长度为m，存放m个指定的插值点处的函数值
'          dz    - Double型一维数组，长度为m，存放m个指定的插值点处的一阶导数值
'          ddz    - Double型一维数组，长度为m，存放m个指定的插值点处的二阶导数值
'  返回值：Double型，指定函数的x(1)到x(n)的定积分值
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INSpline1(n As Integer, x() As Double, y() As Double, dy() As Double, ddy() As Double, m As Integer, t() As Double, z() As Double, dz() As Double, ddz() As Double) As Double
    Dim i As Integer, j As Integer
    Dim h0 As Double, h1 As Double, alpha As Double, beta As Double, g As Double
    ReDim s(n) As Double
    
    ' 初值
    s(1) = dy(1)
    dy(1) = 0#
    h0 = x(2) - x(1)
    
    For j = 2 To n - 1
        h1 = x(j + 1) - x(j)
        alpha = h0 / (h0 + h1)
        beta = (1# - alpha) * (y(j) - y(j - 1)) / h0
        beta = 3# * (beta + alpha * (y(j + 1) - y(j)) / h1)
        dy(j) = -alpha / (2# + (1# - alpha) * dy(j - 1))
        s(j) = (beta - (1# - alpha) * s(j - 1))
        s(j) = s(j) / (2# + (1# - alpha) * dy(j - 1))
        h0 = h1
    Next j
    
    ' 一阶导数值
    For j = n - 1 To 1 Step -1
      dy(j) = dy(j) * dy(j + 1) + s(j)
    Next j
    
    For j = 1 To n - 1
        s(j) = x(j + 1) - x(j)
    Next j
        
    ' 二阶导数值
    For j = 1 To n - 1
        h1 = s(j) * s(j)
        ddy(j) = 6# * (y(j + 1) - y(j)) / h1 - 2# * (2# * dy(j) + dy(j + 1)) / s(j)
    Next j
    
    h1 = s(n - 1) * s(n - 1)
    ddy(n) = 6# * (y(n - 1) - y(n)) / h1 + 2# * (2# * dy(n) + dy(n - 1)) / s(n - 1)
    g = 0#
    For i = 1 To n - 1
        h1 = 0.5 * s(i) * (y(i) + y(i + 1))
        h1 = h1 - s(i) * s(i) * s(i) * (ddy(i) + ddy(i + 1)) / 24#
        g = g + h1
    Next i
    
    ' 插值
    For j = 1 To m
        If (t(j) >= x(n - 1)) Then
            i = n - 1
        Else
            i = 1
            While (t(j) > x(i + 1))
                i = i + 1
            Wend
        End If
        
        h1 = (x(i + 1) - t(j)) / s(i)
        h0 = h1 * h1
        z(j) = (3# * h0 - 2# * h0 * h1) * y(i)
        z(j) = z(j) + s(i) * (h0 - h0 * h1) * dy(i)
        dz(j) = 6# * (h0 - h1) * y(i) / s(i)
        dz(j) = dz(j) + (3# * h0 - 2# * h1) * dy(i)
        ddz(j) = (6# - 12# * h1) * y(i) / (s(i) * s(i))
        ddz(j) = ddz(j) + (2# - 6# * h1) * dy(i) / s(i)
        h1 = (t(j) - x(i)) / s(i)
        h0 = h1 * h1
        z(j) = z(j) + (3# * h0 - 2# * h0 * h1) * y(i + 1)
        z(j) = z(j) - s(i) * (h0 - h0 * h1) * dy(i + 1)
        dz(j) = dz(j) - 6# * (h0 - h1) * y(i + 1) / s(i)
        dz(j) = dz(j) + (3# * h0 - 2# * h1) * dy(i + 1)
        ddz(j) = ddz(j) + (6# - 12# * h1) * y(i + 1) / (s(i) * s(i))
        ddz(j) = ddz(j) - (2# - 6# * h1) * dy(i + 1) / s(i)
    Next j
    
    ' 返回积分值
    INSpline1 = g

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INSpline2
'  功能：  实现第二种边界条件的三次样条函数插值、微商与积分
'  参数：  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
'          dy    - Double型一维数组，长度为n，返回时，存放n个给定点处的一阶导数值y'(i)，i=1,2,…n
'          ddy  - Double型一维数组，长度为n，调用时，ddy(1)存放给定区间的左端点处的二阶导数值，ddy(n)存放
'                 给定区间的右端点处的二阶导数值。返回时，存放n个给定点处的二阶导数值y''(i)，i=1,2,…n
'          m    - Integer型变量，指定插值点的个数
'          t     - Double型一维数组，长度为m，存放m个指定的插值点的值。要求x(1)<t(j)<x(n), j=1,2,…,m-1
'          z     - Double型一维数组，长度为m，存放m个指定的插值点处的函数值
'          dz    - Double型一维数组，长度为m，存放m个指定的插值点处的一阶导数值
'          ddz    - Double型一维数组，长度为m，存放m个指定的插值点处的二阶导数值
'  返回值：Double型，指定函数的x(1)到x(n)的定积分值
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INSpline2(n As Integer, x() As Double, y() As Double, dy() As Double, ddy() As Double, m As Integer, t() As Double, z() As Double, dz() As Double, ddz() As Double) As Double
    Dim i As Integer, j As Integer
    Dim h0 As Double, h1 As Double, alpha As Double, beta As Double, g As Double
    ReDim s(n) As Double
    
    ' 初值
    dy(1) = -0.5
    h0 = x(2) - x(1)
    s(1) = 3# * (y(2) - y(1)) / (2# * h0) - ddy(1) * h0 / 4#
    
    For j = 2 To n - 1
        h1 = x(j + 1) - x(j)
        alpha = h0 / (h0 + h1)
        beta = (1# - alpha) * (y(j) - y(j - 1)) / h0
        beta = 3# * (beta + alpha * (y(j + 1) - y(j)) / h1)
        dy(j) = -alpha / (2# + (1# - alpha) * dy(j - 1))
        s(j) = (beta - (1# - alpha) * s(j - 1))