interpmodule.bas

科学与工程数值计算算法（Visual Basic版） 附赠的光盘包含了本书中全部的源代码

Attribute VB_Name = "InterpModule"
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  功能：  插值算法
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Option Explicit

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INLagrn
'  功能：  用拉格朗日插值公式进行一元全区间不等距插值
'  参数：  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
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INLagrn(n As Integer, x() As Double, y() As Double, t As Double) As Double
    Dim i As Integer, j As Integer, k As Integer, m As Integer
    Dim z As Double, s As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INLagrn = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INLagrn = z
        Exit Function
    End If
    
    If (n = 2) Then
        z = (y(1) * (t - x(2)) - y(2) * (t - x(1))) / (x(1) - x(2))
        INLagrn = z
        Exit Function
    End If
    
    ' 开始插值
    i = 0
    While ((x(i) < t) And (i <= n))
        i = i + 1
    Wend
    
    k = i - 4
    If (k < 0) Then k = 0
    
    m = i + 3
    If (m > n - 1) Then m = n - 1
    
    For i = k To m
        s = 1#
        For j = k To m
            If (j <> i) Then s = s * (t - x(j + 1)) / (x(i + 1) - x(j + 1)) ' 拉格朗日插值公式

        Next j
        z = z + s * y(i + 1)
    Next i
        
    ' 返回结果
    INLagrn = z

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdLagrn
'  功能：  用拉格朗日插值公式进行一元全区间等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          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 INEdLagrn(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
    Dim z As Double, s As Double, xi As Double, xj As Double
    Dim p As Double, q As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INEdLagrn = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INEdLagrn = z
        Exit Function
    End If
        
    If (n = 2) Then
        z = (y(2) * (t - x0) - y(1) * (t - x0 - h)) / h
        INEdLagrn = z
        Exit Function
    End If
    
    ' 开始插值
    If (t > x0) Then
        p = (t - x0) / h
        i = Int(p)
        q = i
        If (p > q) Then i = i + 1
    Else
        i = 0
    End If
    
    k = i - 4
    If (k < 0) Then k = 0
    
    m = i + 3
    If (m > n - 1) Then m = n - 1
    
    For i = k To m
        s = 1#
        xi = x0 + i * h
        For j = k To m
            If (j <> i) Then
                xj = x0 + j * h
                ' 拉格朗日插值公式
                s = s * (t - xj) / (xi - xj)
            End If
        Next j
        z = z + s * y(i + 1)
    Next i
        
    ' 返回结果
    INEdLagrn = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INLagrn3
'  功能：  进行一元三点不等距插值
'  参数：  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
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INLagrn3(n As Integer, x() As Double, y() As Double, t As Double) As Double
    Dim i As Integer, j As Integer, k As Integer, m As Integer
    Dim z As Double, s As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INLagrn3 = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INLagrn3 = z
        Exit Function
    End If
    
    If (n = 2) Then
        z = (y(1) * (t - x(2)) - y(2) * (t - x(1))) / (x(1) - x(2))
        INLagrn3 = z
        Exit Function
    End If
    
    ' 开始插值
    If (t <= x(2)) Then
        k = 0
        m = 2
    Else
        If (t >= x(n - 1)) Then
            k = n - 3
            m = n - 1
        Else
            k = 1
            m = n
            
            While (m - k <> 1)
                i = (k + m) / 2
                If (t < x(i)) Then
                    m = i
                Else
                    k = i
                End If
            Wend
            
            k = k - 1
            m = m - 1
            If (Abs(t - x(k + 1)) < Abs(t - x(m + 1))) Then
                k = k - 1
            Else
                m = m + 1
            End If
        End If
    End If
    
    z = 0#
    For i = k To m
        s = 1#
        For j = k To m
            If (j <> i) Then
                ' 抛物线插值公式
                s = s * (t - x(j + 1)) / (x(i + 1) - x(j + 1))
            End If
        Next j
        
        z = z + s * y(i + 1)
    Next i
    
    ' 返回结果
    INLagrn3 = z
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdLagrn3
'  功能：  进行一元三点等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          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 INEdLagrn3(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
    Dim z As Double, s As Double, xi As Double, xj As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INEdLagrn3 = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INEdLagrn3 = z
        Exit Function
    End If
        
    If (n = 2) Then
        z = (y(2) * (t - x0) - y(1) * (t - x0 - h)) / h
        INEdLagrn3 = z
        Exit Function
    End If

    ' 开始插值
    If (t <= x0 + h) Then
        k = 0
        m = 2
    Else
        If (t >= x0 + (n - 3) * h) Then
            k = n - 3
            m = n - 1
        Else
            i = Int((t - x0) / h) + 1
            If (Abs(t - x0 - i * h) >= Abs(t - x0 - (i - 1) * h)) Then
                k = i - 2
                m = i
            Else
                k = i - 1
                m = i + 1
            End If
        End If
    End If
        
    z = 0#
    For i = k To m
        s = 1#
        xi = x0 + i * h
        
        For j = k To m
            If (j <> i) Then
                xj = x0 + j * h
                ' 抛物线插值公式
                s = s * (t - xj) / (xi - xj)
            End If
        Next j
        
        z = z + s * y(i + 1)
    Next i
    
    ' 返回结果
    INEdLagrn3 = z

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INPq
'  功能：  进行连分式不等距插值
'  参数：  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
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INPq(n As Integer, x() 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, h As Double, b(8) As Double
    
    ' 初值
    z = 0#
    
    ' 特例处理
    If (n < 1) Then
        INPq = z
        Exit Function
    End If
    
    If (n = 1) Then
        z = y(1)
        INPq = z
        Exit Function
    End If
    
    ' 开始插值
    If (n <= 8) Then
        k = 0
        m = n
    Else
        If (t < x(5)) Then
            k = 0
            m = 8
        Else
            If (t > x(n - 4)) Then
                k = n - 8
                m = 8
            Else
                k = 1
                j = n
                While (j - k <> 1)
                    i = (k + j) / 2
                    If (t < x(i)) Then
                        j = i
                    Else
                        k = i
                    End If
                Wend
                
                k = k - 4
                m = 8
            End If
        End If
    End If
    
    b(1) = y(k + 1)
    
    For i = 2 To m
        h = y(i + k)
        l = 0
        j = 1
        
        While ((l = 0) And (j <= i - 1))
            If (Abs(h - b(j)) + 1# = 1#) Then
                l = 1
            Else
                h = (x(i + k) - x(j + k)) / (h - b(j))
            End If
              
            j = j + 1
        Wend
        
        b(i) = h
        
        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 - x(i + k)) / z
    Next i
    
    ' 返回结果
    INPq = z

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INEdPq
'  功能：  进行连分式等距插值
'  参数：  n     - Integer型变量，给定结点的点数