二维插值法,想用P、T联合查W值,请高手指导一下。
以下是数据库内容:P T W
0.11 10 42.05
0.12 20 42.063
0.13 30 42.076
0.14 40 42.089
0.15 50 42.101
0.16 60 42.114
0.17 70 42.127
0.18 80 42.14
0.19 90 42.15
想用P、T查W值,请高手指导一下。
使用插值:
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名: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
form代码如下:
Option Explicit
Dim cn As New ADODB.Connection
Dim rs As New ADODB.Recordset
Private Sub test_Click()
Dim x() As Double, y() As Double, z() As Double ', x As Single, y As Single, z As Single
Dim m As Integer, n As Integer, s As Integer
Dim i As Integer, j As Integer
Dim u As Double, v As Double, w As Double
Dim str As String
cn.Open
rs.Open "H", cn, adOpenKeyset, adLockPessimistic
rs.MoveFirst
On Error Resume Next
For i = 1 To rs.RecordCount Step 1
ReDim Preserve x(i)
x(i) = rs.Fields("P") '取列的数据
rs.MoveNext
'Debug.Print x(i)
Next i
For j = 1 To rs.RecordCount Step 1
ReDim Preserve y(j)
y(j) = rs.Fields("T") '取列的数据
rs.MoveNext
Debug.Print y(j)
Next j
rs.MoveFirst
i = rs.RecordCount
j = rs.Fields.Count - 1
ReDim z1(i, j)
For i = 1 To rs.RecordCount Step 1
j = 1
z1(i, j) = CDbl(rs.Fields(2))
For s = 3 To rs.Fields.Count - 1 '使用循环获取列的信息
j = j + 1
If Not IsNull(rs.Fields(s)) Then z1(i, j) = CDbl(rs.Fields(j + 1))
'Debug.Print z1(i, j)
Next s
rs.MoveNext
Next i
m = rs.RecordCount
n = rs.Fields.Count - 1
cn.Close
u = 0.11
v = 10
' 插值
w = INLagrn2(n, m, x, y, z, u, v)
str = str & "z(" & u & ", " & v & ") = " & w
MsgBox str
End Sub
Private Sub Form_Load()
cn.ConnectionString = "Provider=Microsoft.Jet.OLEDB.4.0;Data Source=" & App.Path & "\main.mdb;Persist Security Info=False"
End Sub