/*************************************************** * Program to demonstrate Newton interpolation * * of Function SIN(X) in double precision * * ------------------------------------------------ * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* * * * C++ Version by J.-P. Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------ * * SAMPLE RUN: * * * * Newton interpolation of SIN(X): * * * * Input X: 0.5 * * * * SIN(X) = 0.47961461 (exact value: 0.47942554) * * Error estimate: -0.00000028 * * * ***************************************************/ #include #include //Label: 100 #define SIZE 15 double X[SIZE],Y[SIZE],Y1[SIZE],Y2[SIZE],Y3[SIZE]; double e,xx,yy; int iv, n; /******************************************************* * Newton divided differences interpolation subroutine * * ---------------------------------------------------- * * Calculates cubic interpolation for a given table. * * iv is the total number of table values. X(i), Y(i) * * are the coordinate table values, Y(i) being the * * dependant variable. The X(i) may be arbitrarily spa- * * ced. xx is the interpolation point which is assumed * * to be in the interval with at least one table value * * to the left, and 3 to the right. If this is violated,* * n will be set to zero. It is assumed that the table * * values are in ascending X(i) order. e is the error * * estimate. The returned value is yy. * *******************************************************/ void Interpol_Newton() { //Labels: e100,e200,e300 int i; double a,b,c; // Check to see if interpolation point is correct n=1; if (xx>=X[1]) goto e100; n=0 ; return; e100: if (xx<=X[iv-3]) goto e200; n=0 ; return; // Generate divided differences e200: for (i=1 ; iX[i]) goto e300; i--; // Begin interpolation a=xx-X[i]; b=a*(xx-X[i+1]); c=b*(xx-X[i+2]); yy=Y[i]+a*Y1[i]+b*Y2[i]+c*Y3[i]; //Calculate next term in the expansion for an error estimate e=c*(xx-X[i+3])*yy/24.0; } void main() { iv=14; printf("\n Newton interpolation of SIN(X):\n"); /* Input sine table ----------------------------------------------------------------- Sine table values from Handbook of mathematical functions by M. Abramowitz and I.A. Stegun, NBS, june 1964 -----------------------------------------------------------------*/ X[1]=0.000; Y[1]=0.00000000; X[2]=0.125; Y[2]=0.12467473; X[3]=0.217; Y[3]=0.21530095; X[4]=0.299; Y[4]=0.29456472; X[5]=0.376; Y[5]=0.36720285; X[6]=0.450; Y[6]=0.43496553; X[7]=0.520; Y[7]=0.49688014; X[8]=0.589; Y[8]=0.55552980; X[9]=0.656; Y[9]=0.60995199; X[10]=0.721; Y[10]=0.66013615; X[11]=0.7853981634; Y[11]=0.7071067812; X[12]=0.849; Y[12]=0.75062005; X[13]=0.911; Y[13]=0.79011709; X[14]=0.972; Y[14]=0.82601466; /*----------------------------------------------------------------- Input interpolation point */ e100: printf("\n"); printf(" Input X: "); scanf("%lf",&xx); Interpol_Newton(); if (n!=0) { printf("\n SIN(X) = %10.8f (exact value: %10.8f)\n",yy,sin(xx)); printf(" Error estimate: %10.8f\n",e); goto e100; } printf("\n"); } // End of file Newton.cpp