常用水文预报算法和计算程序VB版

Public Sub GETtnk(R() As Double, q() As Double, DT As Integer, N As Double, k As Double)

'根据入流和出流时间序列，用矩法算出本场洪水的纳须单位线参数 n 和 k

'R 净雨过程
'Q R相应的出流过程
'DT 计算时段长

Dim MI1 As Double, MI2 As Double, MO1 As Double, MO2 As Double, NI2 As Double, NO2 As Double

Dim X1 As Double, X2 As Double, X3 As Double
Dim i As Integer, XM As Double

Dim NR As Integer, NQ As Integer

NR = UBound(R())
NQ = UBound(q)

i = 0
X1 = 0
X2 = 0
X3 = 0

For i = 1 To NR

X1 = X1 + R(i) * (2 * i - 1)
X2 = X2 + R(i)
X3 = X3 + R(i) * (2 * i - 1) ^ 2

Next
MI1 = X1 / X2 * DT / 2
MI2 = X3 / X2 * DT * DT / 4
X1 = 0
X2 = 0
X3 = 0
For i = 1 To NQ - 1
XM = (q(i) + q(i + 1)) / 2
X1 = X1 + XM * (2 * i - 1)
X2 = X2 + XM
X3 = X3 + XM * (2 * i - 1) ^ 2

Next

MO1 = X1 / X2 * DT / 2
MO2 = X3 / X2 * DT * DT / 4
NI2 = MI2 - MI1 ^ 2
NO2 = MO2 - MO1 ^ 2
N = (MO1 - MI1) ^ 2 / (NO2 - NI2)
k = (NO2 - NI2) / (MO1 - MI1)

End Sub

Public Sub CUNIT(N As Double, KK As Double, DT As Integer, UH() As Double, MAX As Integer, M As Integer)
'//时段单位线计算

Dim IT As Double, TT As Double, ni As Double
Dim JC As Double
Dim X1 As Double, X2 As Double, UNIT As Double, WC As Double

N = Int(N)
UH(1) = 0
IT = 2
TT = DT
10 X1 = 0
X2 = 0

For i = 1 To N
ni = N - 1
JC = 1
If ni = 0 Then GoTo 30
For j = 1 To ni
JC = JC * j
Next
' End If
30 X1 = ((TT + DT) / KK) ^ ni / JC + X1
X2 = (TT / KK) ^ ni / JC + X2

Next
X1 = X1 * Exp(DT / KK)
UH(IT) = Exp(-TT / KK) * (X1 - X2)
If ((UH(IT) < UH(IT - 1) And UH(IT) < 0.0005)) Then
M = IT - 1
UNIT = 0
For i = 1 To M
UNIT = UNIT - UH(i)
Next
WC = UNIT - 1
If Abs(WC) >= 0.0005 Then
UH(M) = UH(M) - WC
If (UH(M) < 0.0005) Then M = M - 1
End If

Else
TT = TT + DT
IT = IT + 1
GoTo 10
End If

End Sub

' 模块名：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
' 函数名：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
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
' 函数名：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