interpmodule.bas

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

        s(j) = s(j) / (2# + (1# - alpha) * dy(j - 1))
        h0 = h1
    Next j
    
    ' 一阶导数值
    dy(n) = (3# * (y(n) - y(n - 1)) / h1 + ddy(n) * h1 / 2# - s(n - 1)) / (2# + dy(n - 1))
    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)) Then
            i = n - 2
        Else
            i = 0
            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
    
    ' 返回积分值
    INSpline2 = g

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INSpline3
'  功能：  实现第三种边界条件的三次样条函数插值、微商与积分
'  参数：  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
'                  要求y(1)=y(n)
'          dy    - Double型一维数组，长度为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 INSpline3(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 y0 As Double, y1 As Double, h0 As Double, h1 As Double, alpha As Double, beta As Double, u As Double, g As Double
    ReDim s(n) As Double
    
    ' 初值
    h0 = x(n) - x(n - 1)
    y0 = y(n) - y(n - 1)
    dy(1) = 0#
    ddy(1) = 0#
    ddy(n) = 0#
    s(1) = 1#
    s(n) = 1#
    
    For j = 2 To n
        h1 = h0
        y1 = y0
        h0 = x(j) - x(j - 1)
        y0 = y(j) - y(j - 1)
        alpha = h1 / (h1 + h0)
        beta = 3# * ((1# - alpha) * y1 / h1 + alpha * y0 / h0)
        If (j < n) Then
            u = 2# + (1# - alpha) * dy(j - 1)
            dy(j) = -alpha / u
            s(j) = (alpha - 1#) * s(j - 1) / u
            ddy(j) = (beta - (1# - alpha) * ddy(j - 1)) / u
        End If
    Next j
    
    ' 二阶导数值
    For j = n - 1 To 2 Step -1
        s(j) = dy(j) * s(j + 1) + s(j)
        ddy(j) = dy(j) * ddy(j + 1) + ddy(j)
    Next j
    
    ' 一阶导数值
    dy(n - 1) = (beta - alpha * ddy(2) - (1# - alpha) * ddy(n - 1)) / (alpha * s(2) + (1# - alpha) * s(n - 1) + 2#)
    For j = 3 To n
        dy(j - 2) = s(j - 1) * dy(n - 1) + ddy(j - 1)
    Next j
    
    dy(n) = dy(1)
    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
        h0 = t(j)
        While (h0 >= x(n))
            h0 = h0 - (x(n) - x(1))
        Wend
        
        While (h0 < x(1))
            h0 = h0 + (x(n) - x(1))
        Wend
        
        i = 1
        While (h0 > x(i + 1))
            i = i + 1
        Wend
        
        u = h0
        h1 = (x(i + 1) - u) / 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 = (u - 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
    
    ' 返回积分值
    INSpline3 = g

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INTqip
'  功能：  实现二元三点插值
'  参数：  n     - Integer型变量，x方向上给定结点的点数
'          m     - Integer型变量，y方向上给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)，要求x(1)<x(2)<...<x(n)
'          y     - Double型一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)，要求y(1)<y(2)<...<y(m)
'          z     - Double型n x m二维数组，存放给定的n x m个结点的函数值z(i,j)，z(i,j) = f(x(i), y(j)), i=1,2,...,n, j=1,2,...,m
'          u     - Double型变量，存放插值点x坐标
'          v     - Double型变量，存放插值点y坐标
'  返回值：Double型，指定函数值f(u, v)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INTqip(n As Integer, m As Integer, x() As Double, y() As Double, z() As Double, u As Double, v As Double) As Double
    Dim nn As Integer, mm As Integer, ip As Integer, iq As Integer, i As Integer, j As Integer, k As Integer, l As Integer
    Dim b(3) As Double, h As Double, w As Double
    
    nn = 3
    If (n <= 3) Then
        ip = 0
        nn = n
    Else
        If (u <= x(2)) Then
            ip = 0
        Else
            If (u >= x(n - 1)) Then
                ip = n - 3
            Else
                i = 1
                j = n
                While (((i - j) <> 1) And ((i - j) <> -1))
                    l = (i + j) / 2
                    If (u < x(l)) Then
                        j = l
                    Else
                        i = l
                    End If
                Wend
                If (Abs(u - x(i)) < Abs(u - x(j))) Then
                    ip = i - 2
                Else
                    ip = i - 1
                End If
            End If
        End If
    End If
    
    mm = 3
    If (m <= 3) Then
        iq = 0
        mm = m
    Else
        If (v <= y(2)) Then
            iq = 0
        Else
            If (v >= y(m - 1)) Then
                iq = m - 3
            Else
                i = 1
                j = m
                While (((i - j) <> 1) And ((i - j) <> -1))
                    l = (i + j) / 2
                    If (v < y(l)) Then
                        j = l
                    Else
                        i = l
                    End If
                Wend
                If (Abs(v - y(i)) < Abs(v - y(j))) Then
                    iq = i - 2
                Else
                    iq = i - 1
                End If
            End If
        End If
    End If
    
    For i = 1 To nn
        b(i) = 0#
        For j = 1 To mm
            h = z(ip + i, iq + j)
            For k = 1 To mm
                If (k <> j) Then
                    ' 二元拉格朗日公式
                    h = h * (v - y(iq + k)) / (y(iq + j) - y(iq + k))
                End If
            Next k
            b(i) = b(i) + h
        Next j
    Next i
    
    w = 0#
    For i = 1 To nn
        h = b(i)
        For j = 1 To nn
            If (j <> i) Then
                ' 二元拉格朗日公式
                h = h * (u - x(ip + j)) / (x(ip + i) - x(ip + j))
            End If
        Next j
        w = w + h
    Next i
    
    ' 返回函数值
    INTqip = w

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：InterpModule.bas
'  函数名：INLagrn2
'  功能：  实现二元全区间插值
'  参数：  n     - Integer型变量，x方向上给定结点的点数
'          m     - Integer型变量，y方向上给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定n x m 个结点x方向上的n个值x(i)，要求x(1)<x(2)<...<x(n)
'          y     - Double型一维数组，长度为m，存放给定n x m 个结点y方向上的m个值y(i)，要求y(1)<y(2)<...<y(m)
'          z     - Double型n x m二维数组，存放给定的n x m个结点的函数值z(i,j)，z(i,j) = f(x(i), y(j)), i=1,2,...,n, j=1,2,...,m
'          u     - Double型变量，存放插值点x坐标
'          v     - Double型变量，存放插值点y坐标
'  返回值：Double型，指定函数值f(u, v)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function INLagrn2(n As Integer, m As Integer, x() As Double, y() As Double, z() As Double, u As Double, v As Double) As Double
    Dim ip As Integer, ipp As Integer, iq As Integer, iqq As Integer, i As Integer, j As Integer, k As Integer, l As Integer
    Dim b(10) As Double, h As Double, w As Double
    
    If (u <= x(1)) Then
        ip = 1
        ipp = 4
    Else
        If (u >= x(n)) Then
            ip = n - 3
            ipp = n
        Else
            i = 1
            j = n
            While (((i - j) <> 1) And ((i - j) <> -1))
                l = (i + j) / 2
                If (u < x(l)) Then
                    j = l
                Else
                    i = l
                End If
            Wend
            ip = i - 3
            ipp = i + 4
        End If
    End If
    
    If (ip < 1) Then
        ip = 1
    End If
    
    If (ipp > n) Then
        ipp = n
    End If
    
    If (v <= y(1)) Then
        iq = 1
        iqq = 4
    Else
        If (v >= y(m)) Then
            iq = m - 3
            iqq = m
        Else
            i = 1
            j = m
            While (((i - j) <> 1) And ((i - j) <> -1))
                l = (i + j) / 2
                If (v < y(l)) Then
                    j = l
                Else
                    i = l
                End If
            Wend
            iq = i - 3
            iqq = i + 4
        End If
    End If
    
    If (iq < 1) Then
        iq = 1
    End If
    
    If (iqq > m) Then
        iqq = m
    End If
    
    For i = ip To ipp
        b(i - ip + 1) = 0#
        For j = iq To iqq
            h = z(i, j)
            For k = iq To iqq
                ' 二元拉格朗日公式
                If (k <> j) Then h = h * (v - y(k)) / (y(j) - y(k))
            Next k
            b(i - ip + 1) = b(i - ip + 1) + h
        Next j
    Next i
    
    w = 0#
    For i = ip To ipp
        h = b(i - ip + 1)
        For j = ip To ipp
          ' 二元拉格朗日公式
          If (j <> i) Then h = h * (u - x(j)) / (x(i) - x(j))
        Next j
        w = w + h
    Next i
    
    ' 返回函数值
    INLagrn2 = w

End Function