'**************************************************** '* 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].* '* * '* ------------------------------------------------ * '* SAMPLE RUN: * '* * '* Newton interpolation of SIN(X): * '* * '* Input X: 0.5 * '* * '* SIN(X) = 0.47961461 (exact value: 0.47942554) * '* Error estimate: -0.00000028 * '* * '**************************************************** defint i-n defdbl a-h,o-z iv=14 dim X(iv+1),Y(iv+1),Y1(iv),Y2(iv-1),Y3(iv-2) cls print print " Newton interpolation of SIN(X):" 'Input table for i=1 to iv read X(i), Y(i) next i '------------------------------------------------------------ 'Sine table values from Handbook of mathematical 'functions by M. Abramowitz and I.A. Stegun, 'NBS, june 1964 data 0.000#,0.00000000#,0.125#,0.12467473# data 0.217#,0.21530095#,0.299#,0.29456472#,0.376#,0.36720285# data 0.450#,0.43496553#,0.520#,0.49688014#,0.589#,0.55552980# data 0.656#,0.60995199#,0.721#,0.66013615# data 0.7853981634#,0.7071067812# data 0.849#,0.75062005#,0.911#,0.79011709# data 0.972#,0.82601466# '------------------------------------------------------------ 'Input interpolation point 100 print input " Input X: ",xx gosub 1000 if n<>0 then print print using " SIN(X) = #.########"; yy; print using " (exact value: #.########)"; SIN(xx) print using " Error estimate: ##.########"; e goto 100 end if print end '******************************************************** '* 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. * '******************************************************** 1000 'Check to see if interpolation point is correct n=1 if xx>=X(1) then goto 1100 n=0 : return 1100 if xx<=X(iv-3) then goto 1200 n=0 : return 'Generate divided differences 1200 for i=1 to iv-1 Y1(i)=(Y(i+1)-Y(i))/(X(i+1)-X(i)) next i for i=1 to iv-2 Y2(i)=(Y1(i+1)-Y1(i))/(X(i+1)-X(i)) next i for i=1 to iv-3 Y3(i)=(Y2(i+1)-Y2(i))/(X(i+1)-X(i)) next i 'Find relevant table interval i=0 1300 i=i+1 if xx>X(i) then goto 1300 i=i-1 '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# return 'End of file Newton.bas