!********************************************************* !* 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.329402327209023 * !* Estimated Error = -1.327032626597764E-010 * !* 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_RATINT 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 RATINT(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 RATINT(XA,YA,N,X,Y,DY) !***************************************************** !* Interpolation or Extrapolation of a Discrete * !* Function By a Quotient of Polynomials * !* ------------------------------------------------- * !* 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,TINY=1.E-25) REAL*8 XA(N),YA(N), X,Y,DY REAL*8 C(NMAX),D(NMAX) REAL*8 DD,H,HH,T,W NS=1 HH=DABS(X-XA(1)) DO I=1,N H=DABS(X-XA(I)) IF(H.EQ.0.D0) THEN Y=YA(I) DY=0.D0 RETURN ELSE IF (H.LT.HH) THEN NS=I HH=H ENDIF C(I)=YA(I) D(I)=YA(I)+TINY !The TINY part is needed to prevent a rare END DO !zero-over-zero condition. Y=YA(NS) NS=NS-1 DO M=1,N-1 DO I=1,N-M W=C(I+1)-D(I) H=XA(I+M)-X !H<>0 (tested before) T=(XA(I)-X)*D(I)/H DD=T-C(I+1) IF(DD.EQ.0.) Pause ' Pole at requested value of X.' DD=W/DD D(I)=C(I+1)*DD C(I)=T*DD END DO IF(2*NS.LT.N-M) THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY END DO RETURN END !end of file tratint.f90