/******************************************************** * 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.25500000e+000 * * Estimated Y value = 3.29402327e-001 * * Estimated Error = -1.32471020e-010 * * Exact Y value = 3.29402327e-001 * * * * ----------------------------------------------------- * * 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]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ********************************************************/ #include #include #define NMAX 25 #define TINY 1e-15 double X[NMAX], Y[NMAX]; double XX, YY, XERR; int N; double FCT(double X) { return (sin(X) - 2.0*cos(X)); } void RATINT(double *XA, double *YA, int N, double X, double *Y, double *DY) { /**************************************************** * Interpolation or Extrapolation of a Discrete * * Function By a Quotient of Polynomials * * ------------------------------------------------- * * INPUTS: * * XA: Table of abscissas (N) * * YA: Table of ordinates (N) * * N: Number of points * * X: Interpolation abscissa value * * OUTPUTS: * * Y: Returned estimation of function for X * * DY: Estimated error for Y * ****************************************************/ double C[NMAX], D[NMAX]; double DD,H,HH,T,W; int I,M,NS; NS=1; HH = fabs(X-XA[1]); for (I=1; I<=N; I++) { H = fabs(X-XA[I]); if (H == 0.0) { *Y=YA[I]; *DY=0.0; return; } else if (H < HH) { NS=I; //index of closest table entry HH = H; } C[I]=YA[I]; D[I]=YA[I] + TINY; //TINY is to prevent a zero-over-zero } //condition. *Y=YA[NS]; NS--; for (M=1; M0 (tested before} T=(XA[I]-X)*D[I]/H; DD=T-C[I+1]; if (DD == 0.0) { printf("\n Error: Pole at requested value of X.\n\n"); return; } DD=W/DD; D[I]=C[I+1]*DD; C[I]=T*DD; } if (2*NS < N-M) *DY=C[NS+1]; else { *DY=D[NS]; NS--; } *Y = *Y + *DY; } } void main() { N = 12; //Number of points //define tables X and Y X[1] = 0.0; Y[1]=FCT(X[1]); X[2] = 0.1; Y[2]=FCT(X[2]); X[3] = 0.2; Y[3]=FCT(X[3]); X[4] = 0.3; Y[4]=FCT(X[4]); X[5] = 0.4; Y[5]=FCT(X[5]); X[6] = 0.5; Y[6]=FCT(X[6]); X[7] = 0.6; Y[7]=FCT(X[7]); X[8] = 0.7; Y[8]=FCT(X[8]); X[9] = 0.8; Y[9]=FCT(X[9]); X[10] = 0.9; Y[10]=FCT(X[10]); X[11] = 1.0; Y[11]=FCT(X[11]); X[12] = 1.1; Y[12]=FCT(X[12]); //define interpolation abscissa XX = 1.255; //call interpolation function RATINT(X,Y,N,XX,&YY,&XERR); //print results printf("\n"); printf(" for X = %15.8e\n", XX); printf(" Estimated Y value = %15.8e\n", YY); printf(" Estimated Error = %15.8e\n", XERR); printf(" Exact Y value = %15.8e\n", FCT(XX)); printf("\n"); } // end of file tratint.cpp