!********************************************************* !* 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.25500000000000 * !* Estimated Y value = 0.329402327224849 * !* Estimated Error = -8.258969146682216E-011 * !* Exact Y value = 0.329402327220005 * !* * !* ----------------------------------------------------- * !* 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]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !********************************************************* PROGRAM TEST_POLINT PARAMETER(NMAX=25) REAL*8 X(NMAX), Y(NMAX), XX, YY, XERR, FCT N = 12 !Number of points !define tables X and Y X(1) = 0.d0; Y(1)=FCT(X(1)) X(2) = 0.1d0; Y(2)=FCT(X(2)) X(3) = 0.2d0; Y(3)=FCT(X(3)) X(4) = 0.3d0; Y(4)=FCT(X(4)) X(5) = 0.4d0; Y(5)=FCT(X(5)) X(6) = 0.5d0; Y(6)=FCT(X(6)) X(7) = 0.6d0; Y(7)=FCT(X(7)) X(8) = 0.7d0; Y(8)=FCT(X(8)) X(9) = 0.8d0; Y(9)=FCT(X(9)) X(10) = 0.9d0; Y(10)=FCT(X(10)) X(11) = 1.d0; Y(11)=FCT(X(11)) X(12) = 1.1d0; Y(12)=FCT(X(12)) !define interpolation abscissa XX = 1.255d0 !call interpolation subroutine CALL POLINT(X,Y,N,XX,YY,XERR) !print results print *,' ' print *,' For X =', XX print *,' Estimated Y value =', YY print *,' Estimated Error =', XERR print *,' Exact Y value =', FCT(XX) print *,' ' END FUNCTION FCT(X) REAL*8 X, FCT FCT=DSIN(X) - 2.D0*DCOS(X) RETURN END SUBROUTINE POLINT(XA,YA,N,X,Y,DY) !***************************************************** !* Polynomial Interpolation or Extrapolation * !* of a Discrete Function * !* ------------------------------------------------- * !* INPUTS: * !* XA: Table of abcissas (N) * !* YA: Table of ordinates (N) * !* N: Number of points * !* X: Interpolating abscissa value * !* OUTPUTS: * !* Y: Returned estimation of function for X * !* DY: Estimated error for Y * !***************************************************** PARAMETER(NMAX=25) REAL*8 XA(N),YA(N), X,Y,DY REAL*8 C(NMAX),D(NMAX) REAL*8 DEN,DIF,DIFT,HO,HP,W NS=1 DIF=DABS(X-XA(1)) DO I=1,N DIFT=DABS(X-XA(I)) IF(DIFT.LT.DIF) THEN NS=I !index of closest table entry DIF=DIFT ENDIF C(I)=YA(I) !Initialize the C's and D's D(I)=YA(I) END DO Y=YA(NS) !Initial approximation of Y NS=NS-1 DO M=1,N-1 DO I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.) Pause 'Error: two identical abcissas)' DEN=W/DEN D(I)=HP*DEN !Update the C's and D's C(I)=HO*DEN END DO IF(2*NS.LT.N-M) THEN !After each column in the tableau XA DY=C(NS+1) !is completed, we decide which correction, ELSE !C or D, we want to add to our accumulating DY=D(NS) !value of Y, i.e. which path to take through NS=NS-1 !the tableau, forking up or down. We do this ENDIF !in such a way as to take the most "straight Y=Y+DY !line" route through the tableau to its apex, END DO !updating NS accordingly to keep track of !where we are. This route keeps the partial RETURN !approximations centered (insofar as possible) END !on the target X.The last DY added is thus the !error indication. !end of file tpolint.f90