Polynomial 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
Given arrays A.1,...,A.N and B.1,...,B.N, and given a value V. For polynomial P. of degree N-1 such that AI=A.I; P.AI=B.I, for I=1,...,N, we estimate value Y=P.V

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 separated by blank - Y, and error estimate Dy

 POLINT: procedure expose A. B. parse arg N, V Ns = 1; Dif = ABS(V - A.1) do I = 1 to N   Dift = ABS(V - A.I)   if Dift < Dif     then do; Ns = I; Dif = Dift; end   C.I = B.I; D.I = B.I end Y = B.Ns; Ns = Ns - 1; do M = 1 to N - 1   do I = 1 to N - M     Ho = A.I - V     IpM = I + M; Hp = A.IpM - V     Ip1 = I + 1; W = C.Ip1 - D.I     Den = Ho - Hp     if Den = 0 then       call ERROR "POLINT: Error -",         "two input A's are",         "(to within roundoff) identical"     Den = W / Den; D.I = Hp * Den     C.I = Ho * Den   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

EXAMPLE
The following program

 N = 5 do J = 1 to 5; A.J = J; end B.1 = 1.8;  B.2 = 2.27; B.3 = 3.18 B.4 = 4.01; B.5 = 4.91 say POLINT(N, 2.25) exit POLINT: procedure expose A. B. ...

displays on the screen 2.48801270 0.0114501953.

CONNECTIONS

Literature
Olehla M., Tiser J. Prakticke pouziti FORTRANu
Nakladatestvi dopravy a spoju Praha 1976
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

 cover contents index main page