'*********************************************************
'* Polynomial Interpolation or Extrapolation *
'* of a Discrete Function F(x) *
'* ----------------------------------------------------- *
'* SAMPLE RUN: *
'* (Example: Function sin(x) - 2*cos(x) is given by 12 *
'* points from x=0 to x=1.1. *
'* Extrapolate for x=1.255). *
'* *
'* For X = 1.255 *
'* Estimated Y value = .3294023272245815 *
'* Estimated Error = -8.273064603451457D-11 *
'* Exact Y value = .3294023272200048 *
'* *
'* ----------------------------------------------------- *
'* REFERENCE: "Numerical Recipes, The Art of Scientific *
'* Computing By W.H. Press, B.P. Flannery, *
'* S.A. Teukolsky and W.T. Vetterling, *
'* Cambridge University Press, 1986" *
'* [BIBLI 08]. *
'* *
'* Basic Release By J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'*********************************************************
'PROGRAM TEST_POLINT
DEFDBL A-H, O-Z
DEFINT I-N
N = 12 'Number of points
DIM X(N), Y(N), C(N), D(N)
'define tables X and Y
X(1) = 0#
X(2) = .1
X(3) = .2
X(4) = .3
X(5) = .4
X(6) = .5
X(7) = .6
X(8) = .7
X(9) = .8
X(10) = .9
X(11) = 1#
X(12) = 1.1
FOR I = 1 TO N
X1 = X(I): GOSUB 500: Y(I) = FCT
NEXT I
' define interpolation abscissa
XX = 1.255#
CLS
' call interpolation procedure
GOSUB 1000 'call POLINT(X,Y,N,XX,YY,DY)
' print results
PRINT
PRINT " For X = "; XX
PRINT " Estimated Y value = "; YY
PRINT " Estimated Error = "; DY
X1 = XX: GOSUB 500
PRINT " Exact Y value = "; FCT
PRINT
END 'of main program
500 'FUNCTION FCT(X1)
FCT = SIN(X1) - 2# * COS(X1)
RETURN
1000 'Subroutine POLINT(X,Y,N,XX,YY,DY)
'*****************************************************
'* Polynomial Interpolation or Extrapolation *
'* of a Discreet Function *
'* ------------------------------------------------- *
'* INPUTS: *
'* X: Table of abscissas (N) *
'* Y: Table of ordinates (N) *
'* N: Number of points *
'* XX: Interpolation abscissa value *
'* OUTPUTS: *
'* YY: Returned estimation of function for X *
'* DY: Estimated error for YY *
'*****************************************************
NS = 1
DIF = ABS(XX - X(1))
FOR I = 1 TO N
DIFT = ABS(XX - X(I))
IF DIFT < DIF THEN
NS = I 'index of closest table entry
DIF = DIFT
END IF
C(I) = Y(I) 'Initialize the C"s and D"s
D(I) = Y(I)
NEXT I
YY = Y(NS) 'Initial approximation of Y
NS = NS - 1
FOR M = 1 TO N - 1
FOR I = 1 TO N - M
HO = X(I) - XX
HP = X(I + M) - XX
W = C(I + 1) - D(I)
DEN = HO - HP
IF DEN = 0# THEN
PRINT
PRINT " Error: two identical abscissas."
STOP 'stop program
END IF
DEN = W / DEN
D(I) = HP * DEN 'Update the C's and D's
C(I) = HO * DEN
NEXT I
IF 2 * NS < N - M THEN 'After each column in the tableau XA is completed,
DY = C(NS + 1) 'we decide which correction, C or D, we want to
ELSE 'add to our accumulating value of Y, i.e. which
DY = D(NS) 'path to take through the tableau, forking up or
NS = NS - 1 'down. We do this in such a way as to take the
END IF 'most "straight line" route through the tableau to
'its apex, updating NS accordingly to keep track
YY = YY + DY 'of where we are. This route keeps the partial
NEXT M 'approximations centered (insofar as possible) on
'the target X.The last DY added is thus the error
'indication.
RETURN
'end of file tpolint.bas