!**************************************************** !* Program to demonstrate Newton interpolation * !* of Function SIN(X) in real*8 precision * !* ------------------------------------------------ * !* Reference: BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* !* * !* F90 Version by J.-P. Moreau, Paris. * !* (www.jpmoreau.fr) * !* ------------------------------------------------ * !* SAMPLE RUN: * !* * !* Newton interpolation of SIN(X): * !* * !* Input X: 0.5 * !* * !* SIN(X) = 0.47961465 (exact value: 0.47942554) * !* Error estimate: -0.00000028 * !* * !**************************************************** PROGRAM Newton !Labels: 50,60,100 parameter(SIZE=14) integer iv, n real*8 X(0:SIZE),Y(0:SIZE) real*8 e,xx,yy iv=14 print *,' ' print *,'Newton interpolation of SIN(X):' ! 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 100 print *,' ' write(*,"(' Input X: ')",advance ='no') read *, xx call Interpol_Newton(iv,n,e,xx,yy,X,Y) if (n.ne.0) then write(*,50) yy, dsin(xx) write(*,60) e goto 100 end if print *,' ' 50 format(/' SIN(X) = ',F10.8,' (exact value: ',F10.8) 60 format(' Error estimate: ',F10.8) stop 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. * !******************************************************** Subroutine Interpol_Newton(iv,n,e,xx,yy,X,Y) !Labels: 100,200,300 parameter(SIZE=14) integer i real*8 a,b,c real*8 X(0:SIZE),Y(0:SIZE),Y1(0:SIZE),Y2(0:SIZE),Y3(0:SIZE) real*8 e,xx,yy ! Check to see if interpolation point is correct n=1 if (xx>=X(1)) goto 100 n=0 ; return 100 if (xx<=X(iv-3)) goto 200 n=0 ; return ! Generate divided differences 200 do i = 1, iv-1 Y1(i)=(Y(i+1)-Y(i))/(X(i+1)-X(i)) end do do i = 1, iv-2 Y2(i)=(Y1(i+1)-Y1(i))/(X(i+1)-X(i)); end do do i = 1, iv-3 Y3(i)=(Y2(i+1)-Y2(i))/(X(i+1)-X(i)); end do ! Find relevant table interval i=0 300 i = i + 1 if (xx>X(i)) goto 300 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.d0 return end ! End of file Newton.f90