?? 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
' 調(diào)用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
' 函數(shù)名：INEdAkima
' 功能： 光滑等距插值
' 參數(shù)： n - Integer型變量，給定結(jié)點的點數(shù)
' h - Integer型變量，等距結(jié)點的步長
' x0 - Double型變量，存放等距n個結(jié)點中第一個結(jié)點的值
' y - Double型一維數(shù)組，長度為n，存放給定的n個等距結(jié)點的函數(shù)值y(i)，y(i) = f(x(i)), i=1,2,...,n
' k - Integer型變量，控制參數(shù)，若k>=0，則只計算第k個子區(qū)間[x(k), x(k+1)]上的三次多項式的系數(shù)
' s1，s2，s3，s4；若k<0，則需要計算指定插值點t處的函數(shù)近似值f(t)，并計算所在子區(qū)間的三
' 次多項式系數(shù)s1，s2，s3，s4
' t - Double型變量，存放指定的插值點的值,若k>=0，此參數(shù)不起作用，可為任意值
' s - Double型一維數(shù)組，長度為5，其中s1，s2，s3，s4返回三次多項式的系數(shù)，s5返回指定插值點t處的
' 函數(shù)近似值f(t)（k<0時）或任意值（k>=0時）
' 返回值：Double型，指定的查指點t的函數(shù)近似值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
' 調(diào)用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
' 函數(shù)名：INSpline1
' 功能： 實現(xiàn)第一種邊界條件的三次樣條函數(shù)插值、微商與積分
' 參數(shù)： n - Integer型變量，給定結(jié)點的點數(shù)
' x - Double型一維數(shù)組，長度為n，存放給定的n個結(jié)點的值x(i)，要求x(1)<x(2)<...<x(n)
' y - Double型一維數(shù)組，長度為n，存放給定的n個等距結(jié)點的函數(shù)值y(i)，y(i) = f(x(i)), i=1,2,...,n
' dy - Double型一維數(shù)組，長度為n，調(diào)用時，dy(1)存放給定區(qū)間的左端點處的一階導(dǎo)數(shù)值，dy(n)存放
' 給定區(qū)間的右端點處的一階導(dǎo)數(shù)值。返回時，存放n個給定點處的一階導(dǎo)數(shù)值y'(i)，i=1,2,…n
' ddy - Double型一維數(shù)組，長度為n，返回時，存放n個給定點處的二階導(dǎo)數(shù)值y''(i)，i=1,2,…n
' m - Integer型變量，指定插值點的個數(shù)
' t - Double型一維數(shù)組，長度為m，存放m個指定的插值點的值。要求x(1)<t(j)<x(n), j=1,2,…,m-1
' z - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的函數(shù)值
' dz - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的一階導(dǎo)數(shù)值
' ddz - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的二階導(dǎo)數(shù)值
' 返回值：Double型，指定函數(shù)的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
' 一階導(dǎo)數(shù)值
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
' 二階導(dǎo)數(shù)值
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
' 函數(shù)名：INSpline2
' 功能： 實現(xiàn)第二種邊界條件的三次樣條函數(shù)插值、微商與積分
' 參數(shù)： n - Integer型變量，給定結(jié)點的點數(shù)
' x - Double型一維數(shù)組，長度為n，存放給定的n個結(jié)點的值x(i)，要求x(1)<x(2)<...<x(n)
' y - Double型一維數(shù)組，長度為n，存放給定的n個等距結(jié)點的函數(shù)值y(i)，y(i) = f(x(i)), i=1,2,...,n
' dy - Double型一維數(shù)組，長度為n，返回時，存放n個給定點處的一階導(dǎo)數(shù)值y'(i)，i=1,2,…n
' ddy - Double型一維數(shù)組，長度為n，調(diào)用時，ddy(1)存放給定區(qū)間的左端點處的二階導(dǎo)數(shù)值，ddy(n)存放
' 給定區(qū)間的右端點處的二階導(dǎo)數(shù)值。返回時，存放n個給定點處的二階導(dǎo)數(shù)值y''(i)，i=1,2,…n
' m - Integer型變量，指定插值點的個數(shù)
' t - Double型一維數(shù)組，長度為m，存放m個指定的插值點的值。要求x(1)<t(j)<x(n), j=1,2,…,m-1
' z - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的函數(shù)值
' dz - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的一階導(dǎo)數(shù)值
' ddz - Double型一維數(shù)組，長度為m，存放m個指定的插值點處的二階導(dǎo)數(shù)值
' 返回值：Double型，指定函數(shù)的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))
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -