/******************************************************** * 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 = -8.25896915e-011 * * 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 double X[NMAX], Y[NMAX]; double XX, YY, XERR; int N; double FCT(double X) { return (sin(X) - 2.0*cos(X)); } void POLINT(double *XA, double *YA, int N, double X, double *Y, double *DY) { /**************************************************** * Polynomial Interpolation or Extrapolation * * of a Discrete Function * * ------------------------------------------------- * * 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 DEN,DIF,DIFT,HO,HP,W; int I,M,NS; NS=1; DIF = fabs(X-XA[1]); for (I=1; I<=N; I++) { DIFT = fabs(X-XA[I]); if (DIFT < DIF) { NS=I; //index of closest table entry DIF = DIFT; } C[I]=YA[I]; //Initialize the C's and D's D[I]=YA[I]; } *Y=YA[NS]; //Initial approximation of Y NS--; for (M=1; M