回复 2楼 xbj_hyml
Public Sub interpolation(Data() As Single, ByVal T As Single, ByRef Y As Single)
Dim TT(10) As Integer
'数据库里边的温度数组
Dim h(9) As Single
Dim u(1 To 9) As Single
Dim λ(1 To 9) As Single
Dim g(10) As Single
Dim i As Integer
Dim j As Integer
Dim yy0 As Single
'边界条件
Dim yy10 As Single
'边界条件
Dim gmm As Single
'中间变量
Dim p As Single
'中间变量
Dim q As Single
'中间变量
Dim m(10) As Single
'解空间
Dim aa(0 To 10, 0 To 10) As Single
'系数矩阵
yy0 = 0
'自然边界条件
yy10 = 0
TT(0) = 20
TT(1) = 100
TT(2) = 200
TT(3) = 300
TT(4) = 400
TT(5) = 500
TT(6) = 600
TT(7) = 700
TT(8) = 800
TT(9) = 900
TT(10) = 1000
For i = 0 To 9 Step 1
h(i) = TT(i + 1) - TT(i)
Next i
For i = 1 To 9 Step 1
u(i) = h(i - 1) / (h(i - 1) + h(i))
Next i
For i = 1 To 9 Step 1
λ(i) = h(i) / (h(i - 1) + h(i))
Next i
For i = 1 To 9
gmm = (Data(i + 1) - Data(i)) / h(i) - (Data(i) - Data(i - 1)) / h(i - 1)
g(i) = 6 * gmm / (h(i) + h(i - 1))
Next i
g(0) = 6 * ((Data(1) - Data(0)) / h(0) - yy0) / h(0)
g(10) = 6 * (yy10 - (Data(10) - Data(9)) / h(9)) / h(9)
For i = 0 To 10
'系数矩阵清零
For j = 0 To 10
aa(i, j) = 0
Next j
Next i
aa(0, 0) = 2
'为系数矩阵赋值
aa(0, 1) = 1
For i = 1 To 9
aa(i, i) = 2
aa(i, i - 1) = u(i)
aa(i, i + 1) = λ(i)
Next i
aa(10, 9) = 1
aa(10, 10) = 2
For i = 1 To 10
'矩阵的对角化
p = aa(i - 1, i - 1) / aa(i, i - 1)
g(i) = g(i) * p - g(i - 1)
For j = 0 To 10
aa(i, j) = aa(i, j) * p - aa(i - 1, j)
Next j
Next i
m(10) = g(10) / aa(10, 10)
For i = 9 To 0 Step -1
q = 0
For j = 10 To i + 1 Step -1
q = q + aa(i, j) * m(j)
Next j
m(i) = (g(i) - q) / aa(i, i)
Next i
For i = 0 To 9
If T >= TT(i) And T <= TT(i + 1) Then
Y = (m(i) * (TT(i + 1) - T) ^ 3) / (6 * h(i)) + (m(i + 1) * (T - TT(i)) ^ 3) / (6 * h(i)) + (Data(i) - (m(i) * h(i) ^ 2) / 6) * (TT(i + 1) - T) / h(i) + (Data(i + 1) - (m(i + 1) * h(i) ^ 2) / 6) * (T - TT(i)) / h(i)
End If
Next i
End Sub
这是interpolation