三次样条插值的代码, 求高手帮助~
求 三次样条插值的 vfp代码~ 如, 求b=67.34759时, 对应的 a插值( 是三次样条插值, 不是拉格朗日或分段等其它插值; 不要用到matlab的函数 )~ 不胜感谢~
a b
=== =========
0.5 18.453112
1 26.708874
1.5 33.719939
2 40.313374
2.5 46.442608
3 52.513034
4 64.200686
5 75.733463
6 87.259722
7 98.840930
8 110.455645
9 122.262716
10 134.314002
11 146.523455
12 159.029723
13 171.733376
14 184.845683
15 198.282892
16 211.876066
17 225.987211
18 240.518311
19 255.159935
20 270.492911
21 285.997688
22 302.161308
附: 网上下载的三次样条插值的 FORTRAN/vb 源代码
三次样条插值的FORTRAN源代码
C SPLINE3(N1,T1,X1,N2,T2,X2,D1,D2,C)
C ROLE : CALCULATE INTERPOLATION,DIFFERENTIAL(ONE DEGREE AND TWO DEGREE),CALCULUS
C METHOD: THREE SPLINE INTERPOLATION
C INPUT : N1,T1(N1),X1(N1)=NUMBER OF X1, TIME INDEX OF X1, ORIGINAL DATA
C INPUT : N2,T2(N2)=NUMBER OF X2, TIME INDEX OF X2
C OUTPUT: X2(N2)=INTERPOLATION DATA AT TIME T2 FROM X1(N1)
C OUTPUT: D1(N2),D2(N2),C(N2)=ONE DEGREE DIFFENTIAL,TWO DEGREE DIFFENTIAL,CALCULUS OF VECTOR X2
C INTEGER: N1,N2
C REAL*8 : T1,X1,T2,X2,D1,D2,C
c n = the number of the data samples
c m = the number of the interpolation points you want to produce
c x(n)= epochs of sample points
c y(n)= value of function at samples points
c t(m)= epochs of the interpolation points
c u(m)= values of function at interpolation points
c uu= the first degree differential(dao suo)
c v = the second degree differential(dao shu)
c sum= calculus(ji fen)
c a(n),b(n),c(n),d(n)= the working array
c reference: 丁月蓉，天文数据处理(for interpolation and differential)； 徐士良，FORTRAN常用算法程序集(for calculus)。
subroutine spline3(n,x,y,M,t,u,uu,V,sum)
implicit double precision (a-h,o-z)
dimension x(n),y(n),t(m),u(m),a(n),b(n),c(n),d(n),uu(m),V(m)
dimension sum(m)
a(1)=0.0d0
d(1)=0.0d0
d(n)=0.0d0
c(n)=0.0d0
a(n)=1.0d0
b(1)=1.0d0
c(1)=-1.0d0
b(n)=-1.0d0
l=n-1
do 5 i=2,l
a(i)=(x(i)-x(i-1))/6.0d0
c(i)=(x(i+1)-x(i))/6.0d0
b(i)=2.0d0*(a(i)+c(I))
5 d(i)=(y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1))
c(1)=c(1)/b(1)
d(1)=d(1)/b(1)
do 10 i=2,n
c(i)=c(i)/(b(i)-a(i)*c(i-1))
10 d(i)=(d(i)-a(i)*d(i-1))/(b(i)-a(I)*c(i-1))
a(n)=d(n)
do 15 i=1,l
j=n-i
15 a(j)=d(j)-c(j)*a(j+1)
do 30 j1=1,m
f=t(j1)
do 20 j2=1,n-1
if(x(j2).le.f.and.f.le.x(j2+1)) goto 25
20 continue
goto 30
25 e=x(j2+1)-x(j2)
rr=(a(j2)*(x(j2+1)-f)**3+a(j2+1)*(f-x(j2))**3)/6.0d0/e
ss=(x(j2+1)-f)*(y(j2)/e-a(j2)*e/6.0d0)
tt=(f-x(j2))*(y(j2+1)/e-a(j2+1)*e/6.0d0)
aa=(a(j2+1)*(f-x(j2))**2)/2.0d0/e
dd=(a(j2)*(x(j2+1)-f)**2)/2.0d0/e
bb=(y(j2+1)-y(j2))/e
cc=(a(j2+1)-a(j2))*e/6.0d0
uu(j1)=aa+bb-cc-dd
u(j1)=rr+ss+tt
V(j1)=(a(j2)*(x(j2+1)-f)+a(j2+1)*(f-x(j2)))/e
30 continue
DO 33 I=1,M
do 33 j3=1,I-1
if(j3.ne.m) then
e=t(j3+1)-t(j3)
else
e=t(m)-t(m-1)
endif
sum(I)=sum(I)+.5d0*e*(u(j3+1)+u(j3))-e**3*(V(j3)+V(j3+1))/24.0d0
33 CONTINUE
CPRINT *, U
return
end
vb三次样条插值函数绘图的源代码(绘图部分不需要, 只需 求插入值的代码即可)
Dim X(1000) As Single, Y(1000) As Single
Dim u1(0 To 80000) As Single, v1(0 To 80000) As Single
Dim num As Long
Dim t As Integer
Private Declare Sub Sleep Lib "kernel32" (ByVal dwMilliseconds As Long)
Dim de As Integer
Dim ToInit As Boolean
Dim DownX As Single, DownY As Single
Sub Drawposi(Index As Integer)
Me.Picture1.ForeColor = 0
Me.Picture1.Line (0, Y(Index))-(1024, Y(Index))
Me.Picture1.Line (X(Index), 0)-(X(Index), 770)
End Sub
Function hypot(ByVal X As Single, ByVal Y As Single)
hypot = Sqr(X ^ 2 + Y ^ 2)
End Function
Sub MovePic(Index As Integer)
Dim i As Integer
X(Index) = Picture2(Index).Left + 4
Y(Index) = Picture2(Index).Top + 4
lblX.Caption = "X: " + CStr(CInt(X(Index)))
lblY.Caption = "Y: " + CStr(CInt(Y(Index)))
lblX.Refresh
lblY.Refresh
Me.Picture1.Cls
Me.Picture1.ForeColor = QBColor(10)
For i = 0 To t - 1
Me.Picture1.CurrentX = X(i) + 4
Me.Picture1.CurrentY = Y(i) + 4
Me.Picture1.Print i
Next i
End Sub
Private Sub Command1_Click()
Dim i As Long
'Picture1.Scale (0, 0)-(640, 550)
DrawWidth = 3
Picture1.Cls
'If Check1.Value Then Command2_Click
'X(0) = 1
'Y(0) = 1
'X(t - 1) = 638
'Y(t - 1) = 548
Picture1.ForeColor = QBColor(10)
For i = 0 To t - 1
Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B
Picture1.Print i
Next i
Picture1.ForeColor = QBColor(12)
DrawWidth = 1
tspLine t - 1, 2, 0, 0, 0, 0
Picture1.PSet (u1(0), v1(0))
For i = 1 To num - 1
Picture1.Line -(u1(i), v1(i))
'For de = 1 To 12000: Next de 'Sleep 1
Next i
Picture1.ForeColor = QBColor(10)
For i = 0 To t - 1
Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B
Picture1.Print i
Next i
End Sub
Private Sub Command2_Click()
Dim i As Integer
Randomize Timer
ToInit = Not ToInit
If ToInit Then
= False
= "结束初始化"
Me.Cls
For i = 1 To t - 1
Load Me.Picture2(i)
Next i
For i = 0 To t - 1
Picture2(i).Left = X(i) - 4
Picture2(i).Top = Y(i) - 4
Picture2(i).Visible = True
Next i
Picture1.Cls
Me.Picture1.ForeColor = QBColor(10)
For i = 0 To t - 1
Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B
Picture1.Print i
Next i
Else
= True
= "开始初始化"
For i = 1 To t - 1
Unload Me.Picture2(i)
Next i
Me.Picture2(0).Visible = False
End If
Exit Sub
For i = 0 To t
X(i) = Rnd(1) * 500 + Rnd(1) * 50 + 12
Y(i) = Rnd(1) * 400 + Rnd(1) * 100 + 12
'X(i) = i * 20 + Rnd(1) * 10 + 12
'Y(i) = i * 10 + Rnd(1) * 300 + 22 - Rnd(1) * 200
Next i
End Sub
Sub tspLine(ByVal n As Integer, ByVal ch As Integer, ByVal tx1 As Single, ByVal tx2 As Single, ByVal ty1 As Single, ByVal ty2 As Single)
Dim a(1000) As Single, b(1000) As Single, c(1000) As Single, dX(1000) As Single, dY(1000) As Single
Dim qx(1000) As Single, qy(1000) As Single
Dim tt As Single, bx3 As Single, bx4 As Single, by3 As Single, by4 As Single
Dim cx As Single, cy As Single, t(1000) As Single, px(1000) As Single, py(1000) As Single
Dim u(3000) As Single, v(3000) As Single, i As Integer
num = 0
For i = 1 To n
t(i) = hypot(X(i) - X(i - 1), Y(i) - Y(i - 1))
Next i
Select Case ch
Case 0 '抛物条件
u(0) = (X(1) - X(0)) / t(1): u(1) = (X(2) - X(1)) / t(2)
u(2) = (u(1) - u(0)) / (t(2) + t(1))
tx1 = u(0) - u(2) * t(1)
u(0) = (Y(1) - Y(0)) / t(1): u(1) = (Y(2) - Y(1)) / t(2)
u(2) = (u(1) - u(0)) / (t(2) + t(1))
ty1 = u(0) - u(2) * t(1)
u(0) = (X(n) - X(n - 1)) / t(n): u(1) = (X(n - 1) - X(n - 2)) / t(n - 1)
u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))
tx2 = u(0) + u(2) * t(n)
u(0) = (Y(n) - Y(n - 1)) / t(n): u(1) = (Y(n - 1) - Y(n - 2)) / t(n - 1)
u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))
ty2 = u(0) + u(2) * t(n)
Case 1 '夹持条件
a(0) = 1: c(0) = 0: dX(0) = tx1: dY(0) = ty1
a(n) = 1: b(n) = 0: dX(n) = tx2: dY(n) = ty2
Case 2 '自由条件
a(0) = 2: c(0) = 1
dX(0) = 3 * (X(1) - X(0)) / t(1): dY(0) = 3 * (Y(1) - Y(0)) / t(1)
a(n) = 2: b(n) = 1
dX(n) = 3 * (X(n) - X(n - 1)) / t(n): dY(n) = 3 * (Y(n) - Y(n - 1)) / t(n)
Case 3 '循环条件
a(0) = 2: c(0) = 1
dX(0) = 3 * (X(1) - X(0)) / t(1) - (t(1) * (X(2) - X(1)) / t(2) - X(1) + X(0)) / (t(1) + t(2))
dY(0) = 3 * (Y(1) - Y(0)) / t(1) - (t(1) * (Y(2) - Y(1)) / t(2) - Y(1) + Y(0)) / (t(1) + t(2))
a(n) = 2: b(n) = 1
dX(n) = 3 * (X(n) - X(n - 1)) / t(n)
dX(n) = dX(n) + (X(n) - X(n - 1) - t(n) * (X(n - 1) - X(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))
dY(n) = 3 * (Y(n) - Y(n - 1)) / t(n)
dY(n) = dY(n) + (Y(n) - Y(n - 1) - t(n) * (Y(n - 1) - Y(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))
End Select
'计算方程组系数阵和常数阵
For i = 1 To n - 1
a(i) = 2 * (t(i) + t(i + 1)): b(i) = t(i + 1): c(i) = t(i)
dX(i) = 3 * (t(i) * (X(i + 1) - X(i)) / t(i + 1) + t(i + 1) * (X(i) - X(i - 1)) / t(i))
dY(i) = 3 * (t(i) * (Y(i + 1) - Y(i)) / t(i + 1) + t(i + 1) * (Y(i) - Y(i - 1)) / t(i))
Next i
'采用追赶法解方程组
c(0) = c(0) / a(0)
For i = 1 To n - 1
a(i) = a(i) - b(i) * c(i - 1): c(i) = c(i) / a(i)
Next i
a(n) = a(n) - b(n) * c(i - 1)
qx(0) = dX(0) / a(0): qy(0) = dY(0) / a(0)
For i = 1 To n
qx(i) = (dX(i) - b(i) * qx(i - 1)) / a(i)
qy(i) = (dY(i) - b(i) * qy(i - 1)) / a(i)
Next i
px(n) = qx(n): py(n) = qy(n)
For i = n - 1 To 0 Step -1
px(i) = qx(i) - c(i) * px(i + 1)
py(i) = qy(i) - c(i) * py(i + 1)
Next i
'计算曲线上点的坐标
For i = 0 To n - 1
bx3 = (3 * (X(i + 1) - X(i)) / t(i + 1) - 2 * px(i) - px(i + 1)) / t(i + 1)
bx4 = ((2 * (X(i) - X(i + 1)) / t(i + 1) + px(i) + px(i + 1)) / t(i + 1)) / t(i + 1)
by3 = (3 * (Y(i + 1) - Y(i)) / t(i + 1) - 2 * py(i) - py(i + 1)) / t(i + 1)
by4 = ((2 * (Y(i) - Y(i + 1)) / t(i + 1) + py(i) + py(i + 1)) / t(i + 1)) / t(i + 1)
tt = 0
While (tt <= t(i + 1))
cx = X(i) + (px(i) + (bx3 + bx4 * tt) * tt) * tt
cy = Y(i) + (py(i) + (by3 + by4 * tt) * tt) * tt
u1(num) = cx: v1(num) = cy: num = num + 1: tt = tt + 0.5
Wend
u1(num) = X(i + 1): v1(num) = Y(i + 1): num = num + 1
Next i
End Sub
Private Sub Form_Load()
Dim i As Integer
t = 30
ToInit = False
'Picture1.Scale (0, 0)-(640, 550)
Randomize Timer
= "开始初始化"
For i = 0 To t
X(i) = Rnd(1) * 500 + Rnd(1) * 50 + 12
Y(i) = Rnd(1) * 400 + Rnd(1) * 100 + 12
Next i
For i = 0 To t
X(i) = i * 30 + 20
Y(i) = i * 20 + 20
Next i
'Me.Picture1.Picture = LoadPicture("c:\my documents\MenuBack.bmp")
Me.Picture1.BackColor = QBColor(0)
End Sub
Private Sub Form_Resize()
On Error Resume Next
Me.Picture1.Height = Me.ScaleHeight - 40
End Sub
Private Sub Form_Unload(Cancel As Integer)
End
End Sub
Private Sub Picture2_MouseDown(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)
On Error Resume Next
If Button = 1 Then
DownX = X
DownY = Y
Picture2(Index).ZOrder 0
Picture2(Index - 1).BackColor = QBColor(12)
Picture2(Index + 1).BackColor = QBColor(12)
lblX.Caption = "X: " + CStr(CInt(Picture2(Index).Left + 4))
lblY.Caption = "Y: " + CStr(CInt(Picture2(Index).Top + 4))
lblX.Refresh
lblY.Refresh
MovePic Index
Drawposi Index
End If
End Sub
Private Sub Picture2_MouseMove(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)
If Button = 1 Then
Picture2(Index).Left = Picture2(Index).Left - DownX + X
Picture2(Index).Top = Picture2(Index).Top - DownY + Y
MovePic Index
Command1_Click
Drawposi Index
End If
End Sub
Private Sub Picture2_MouseUp(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)
On Error Resume Next
If Button = 1 Then
DownX = X
DownY = Y
Picture2(Index - 1).BackColor = QBColor(15)
Picture2(Index + 1).BackColor = QBColor(15)
'MovePic Index
lblX.Caption = "X:"
lblY.Caption = "Y:"
lblX.Refresh
lblY.Refresh
Command1_Click
End If
End Sub
[ 本帖最后由 茵梦湖 于 2010-6-4 17:18 编辑 ]
