Rational function interpolation and extrapolation
PROBLEM
We know the value of a function at a set of points A.1<A.2<...<A.N. Estimate the value of a function for arbitrary V.
ALGORITHM
For the interpolation problem (A.1<=V<A.N), a rational function is constructed so as to go through a chosen set of tabulated function values. J. Stoer and R. Bulirsh have found an algorithm for the extrapolation problem (V is outside range).
IMPLEMENTATION
Unit: internal function
Global variables: ascending sequence of values A.1,...,A.N, array B.
Parameters: positive integer N, arbitrary value V
Returns: Couple of values Y, and error estimate Dy
RATINT: procedure expose A. B.
parse arg N, V
Tiny = 1E-25; Ns = 1; Hh = ABS(V - A.1)
do I = 1 to N
H = ABS(V - A.I)
if H = 0
then do
Y = B.I; Dy = 0
return Y Dy
end
else if H < Hh
then do; Ns = I; Hh = H; end
C.I = B.I; D.I = B.I + Tiny
end
Y = B.Ns; Ns = Ns - 1
do M = 1 to N - 1
do I = 1 to N - M
Ip1 = I + 1; W = C.Ip1 - D.I
IpM = I + M; H = A.IpM - V
T = (A.I - V) * D.I / H
Dd = T - C.Ip1
if Dd = 0 then
call ERROR "RATINT: Error -",
"the interpolating function has",
"a pole at the requested value V"
Dd = W / Dd; D.I = C.Ip1 * Dd; C.I = T * Dd
end
if 2 * Ns < (N - M)
then do; Nsp1 = Ns + 1; Dy = C.Nsp1; end
else do; Dy = D.Ns; Ns = Ns - 1; end
Y = Y + Dy
end
return Y Dy
ERROR: say ARG(1); exit
|
CONNECTIONS
Literature
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992
