'********************************************************* '* 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 = .3294023272088324 * '* Estimated Error = -1.327291621149679D-10 * '* 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 RATINT(X,Y,N,XX,YY,DY) '***************************************************** '* Interpolation or Extrapolation of a Discrete * '* Function By a Quotient of Polynomials * '* ------------------------------------------------- * '* 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 * '***************************************************** TINY = .000000000000001# NS = 1 HH = ABS(XX - X(1)) FOR I = 1 TO N H = ABS(XX - X(I)) IF H = 0# THEN Y = Y(I) DY = 0# STOP ELSEIF H < HH THEN NS = I 'index of closest table entry HH = H END IF C(I) = Y(I) D(I) = Y(I) + TINY '{TINY is to prevent a zero-over-zero NEXT I 'condition. YY = Y(NS) NS = NS - 1 FOR M = 1 TO N - 1 FOR I = 1 TO N - M W = C(I + 1) - D(I) H = X(I + M) - XX 'H<>0 (tested before) T = (X(I) - XX) * D(I) / H DD = T - C(I + 1) IF DD = 0# THEN PRINT PRINT " Error: Pole at requested value of X." STOP 'stop program END IF DD = W / DD D(I) = C(I + 1) * DD C(I) = T * DD NEXT I IF 2 * NS < N - M THEN DY = C(NS + 1) ELSE DY = D(NS) NS = NS - 1 END IF YY = YY + DY NEXT M RETURN 'end of file tratint.bas