Interpolation lets you estimate the value of a process between two known tabulated values. Suppose that values of a process y(t) are tabulated at tn with a constant interval of dt = tn+1 - tn. To find the value of y(tn + u·dt), where
linear interpolation gives0 < u < 1,
Often, this is good enough. Sometimes greater accuracy is needed. Given the tabulated values of a process y(t) at six evenly spaced points t1 ... t6; we want to evaluate y(t) at t3 + u·dt, where dx is the interval between samples. In this situation, we can use nonlinear or polynomial interpolation to obtain a more precise value for y(t). A theorem attributed to Weierstrass, which states that every continuous function in an interval (a,b) can be represented in that interval to any desired accuracy by a polynomial, justifies the use of polynomial interpolation. Several methods work, including:y = y(tn) + u·[y(tn+1) - y(tn)]
These methods all use a table of differences for equal-interval interpolation. The first two columns are tn and y(tn). The entries in the remaining columns are the differences between the two entries in the column immediately to the left. Bessel's finite difference interpolation uses a triangular table with entries as shown below:
t1 y1
dA1.5
t2 y2 dB2
dA2.5 dC2.5
t3 y3 dB3 dD3
dA3.5 dC3.5 dE3.5
t4 y4 dB4 dD4
dA4.5 dC4.5
t5 y5 dB5
dA5.5
t6 y6
For the above table, with six samples, Bessel's formula is:
Wherey = 0.5*( y3 + y4 ) + B1(u)*dA3.5 + B2(u)*(dB3 + dB4) + B3(u)*dC3.5 + B4(u)*(dD3 + dD4) + B5(u)*dE3.5
Bessel's coefficients defined as follows:u = (t - t3) / (t4 - t3)
u - 1/2 (u + n/2 - 3/2)!
Bn(u) = ------- · --------------------------------- n is odd
n (u + n/2 - 3/2 - n + 1)! (n - 1)!
1 (u + n/2 - 1)!
Bn(u) = --- · --------------------- n is even
2 (u + n/2 - 1 - n)! n!
The first few coefficients are:
Interpolation formula taken as far as B1 is identical to linear interpolation. The addition of higher terms fits a curve to more and more points on both sides of the desired value to give a more accurate value for y(t).B1(u) = u - 0.5 B2(u) = 0.5·u·(u - 1)/2 B3(u) = (u - 0.5)·u·(u - 1)/6 B4(u) = 0.5·(u + 1)·u·(u - 1)·(u - 2)/24 B5(u) = (u - 0.5)·(u + 1)·u·(u - 1)·(u - 2)/120 B6(u) = 0.5·(u + 2)·(u + 1)·u·(u - 1)·(u - 2)·(u - 3)/720
Below is a sample program showing an implementation of Bessel interpolation to calculate the position of the moon using positions tabulated daily.
Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.
type table : array[6] of real
program
var ras, dec, dst : table
% Moon : Apparent Geocentric Positions, True Equator and Equinox of Date
% UT1 : h m s h m s o ' " km
% ------------------------------------------------------------------------
% 2004 Apr 17 00:00:00.0 23 57 08.676 - 4 16 29.27 389306.996
% 2004 Apr 18 00:00:00.0 0 42 27.849 + 1 37 35.17 392652.502
% 2004 Apr 19 00:00:00.0 1 27 20.166 + 7 22 35.29 395849.635
% 2004 Apr 20 00:00:00.0 2 12 38.939 + 12 45 44.69 398812.496
% 2004 Apr 21 00:00:00.0 2 59 08.486 + 17 34 58.11 401416.265
% 2004 Apr 22 00:00:00.0 3 47 20.082 + 21 38 41.48 403502.112
ras[1] := ras2real( 23, 57, 08.676 ) - 24 % due to discontinuity below
ras[2] := ras2real( 0, 42, 27.849 )
ras[3] := ras2real( 1, 27, 20.166 )
% result between element 3 and 4
ras[4] := ras2real( 2, 12, 38.939 )
ras[5] := ras2real( 2, 59, 08.486 )
ras[6] := ras2real( 3, 47, 20.082 )
dec[1] := dec2real( - 4, 16, 29.27)
dec[2] := dec2real( + 1, 37, 35.17)
dec[3] := dec2real( + 7, 22, 35.29)
% result between element 3 and 4
dec[4] := dec2real( +12, 45, 44.69)
dec[5] := dec2real( +17, 34, 58.11)
dec[6] := dec2real( +21, 38, 41.48)
dst[1] := 389306.996
dst[2] := 392652.502
dst[3] := 395849.635
% result between element 3 and 4
dst[4] := 398812.496
dst[5] := 401416.265
dst[6] := 403502.112
var y : real
var u : real := time2real( 18, 0, 0 )
y := interpolate( ras, u )
put "RA =", ha2str( y )...
y := interpolate( dec, u )
put ", Dec =", dec2str( y )...
y := interpolate( dst, u )
put ", Dist =", y:12:3, " km"
end program
% interpolate from six element table using Bessel's formula
% u - fraction between elements 3 and 4 of table
% y - table of six values assumed to be spaced evenly and to be precise
function interpolate( var y : table, u : real ) : real
var rv : real
% first level of divided differences
var dA15 : real := y[2] - y[1]
var dA25 : real := y[3] - y[2]
var dA35 : real := y[4] - y[3]
var dA45 : real := y[5] - y[4]
var dA55 : real := y[6] - y[5]
% second level
var dB2 : real := dA25 - dA15
var dB3 : real := dA35 - dA25
var dB4 : real := dA45 - dA35
var dB5 : real := dA55 - dA45
% third level
var dC25 : real := dB3 - dB2
var dC35 : real := dB4 - dB3
var dC45 : real := dB5 - dB4
% fourth level
var dD3 : real := dC35 - dC25
var dD4 : real := dC45 - dC35
% fifth level
var dE35 : real := dD4 - dD3
rv := 0.5*(y[3] + y[4]) + B1(u)*dA35 + B2(u)*(dB3 + dB4)
+ B3(u)*dC35 + B4(u)*(dD3 + dD4) + B5(u)*dE35
return rv
end function
% first bessel coefficient
function B1( u : real ) : real
return (u - 0.5)
end function
% second bessel coefficient
function B2( u : real ) : real
return 0.5*u*(u - 1)/2
end function
% third bessel coefficient
function B3( u : real ) : real
return (u - 0.5)*u*(u - 1)/6
end function
% fourth bessel coefficient
function B4( u : real ) : real
return 0.5*(u + 1)*u*(u - 1)*(u - 2)/24
end function
% fifth bessel coefficient
function B5( u : real ) : real
return (u - 0.5)*(u + 1)*u*(u - 1)*(u - 2)/120
end function
% convert (hr min sec) to fraction of day
function time2real( hr, min, sec : int ) : real
return (hr + min/60 + sec/3600)/24
end function
% convert (hr min sec) to right ascension
function ras2real( hr, min : int, sec : real ) : real
return hr + min/60 + sec/3600
end function
% convert (deg min sec) to declination
function dec2real( deg, min : int, sec : real ) : real
var rv : real
if deg > 0 then
rv := deg + min/60 + sec/3600
else
rv := -( -deg + min/60 + sec/3600 )
end if
return rv
end function
% converts hour angle into hour angle string
function ha2str( ra : real ) : string
assert (0 <= ra) and (ra < 24) % assure valid range
var h : int := floor( ra )
ra := ra - h
var m : int := floor( 60*ra )
ra := ra - m / 60
var s : real := 3600*ra
return intstr( h, 3 ) & "h " & intstr( m, 2 ) & "m " & frealstr( s, 5, 3 ) & "s"
end function
% converts declination angle in degrees into string
function dec2str( x : real ) : string
assert (-90 <= x) and (x <= +90) % assure valid range
var dec : real := abs( x ) % work with absolute value
var d : int := floor( dec )
dec := dec - d
var m : int := floor( 60*dec )
dec := dec - m / 60
var s : real := 3600*dec
return intstr( sgn(x)*d, 3 ) & "° " & intstr( m, 2 ) & "' " & frealstr( s, 4, 2 ) & "\""
end function
% returns absolute value of real number
function abs( x : real ) : real
if x < 0 then
return -x
end if
return x
end function
% returns sign of real number as int
function sgn( x : real ) : int
if x < 0 then
return -1
end if
return +1
end function
RA = 2h 1m 14.150s, Dec = 11° 27' 40.49", Dist = 398099.951 km
AbCd Classics - free on-line books