'********************************************************* '* 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