This is POLYL in for\lagrange1.for
Where
By virtue of the fact that one term is always left out of the product that covers N values in which x always appears once, it is easy to see that Lm is a polynomial of degree N-1 and that sums of the Lm values times constants will also be of degree N-1. The value of Lk(xm¹k) where xm¹k is any of the data points other than the k’th one will have a zero in the product for the j=m term so that Lk(xm¹k)=0. The value of Lk(xm=k) is quite different. Each term in the product is of the form and by construction the dangerous term with xj=xk has been left out so that Lk(xm=k)=1.
The Lagrange interpolation implicitly uses the function values to construct the N-1 derivatives in the region with an error equal to the N the derivative anywhere in the region and then uses these to form the polynomial approximating function. This means that the error term in the Lagrange polynomial approximation for an internal value of x is where x is any point between the first and last point and w is on the order of the size of the region. For values of x outside the region w becomes on the order of the distance between x and the most distant point in the interpolation and x is any point in that same region.
Figure 1 shows point used by PolyL
and definition of mbeg, nbeg
The data shown is as returned by BLI. The x axis contains values of xdat, while the y axis contains values of fdat. The 1, 2, …, 5 refers to the value of j in . In the lagrange code these x values are referred to as J+MBEG. The integration of interest is from x=nbeg to x=nbeg+1.
Where
The value nbeg is returned by locate, when it is called as CALL LOCATE(X,XDAT,NDAT,NBEG).
SUBROUTINE UELAG(NL,X,MBEG,NSKIP,ALAG,XDAT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ALAG(*),XDAT(*)
DO M=1,NL
IF(M+MBEG.EQ.NSKIP)THEN
ç NSKIP is for error
analysis.
ALAG(M)=0
ELSE
ALAG(M)=1
DO J=1,NL
IF(J.NE.M.AND.J+MBEG.NE.NSKIP)THEN
ALAG(M)=ALAG(M)*(X-XDAT(J+MBEG))/(XDAT(M+MBEG)-
2
XDAT(J+MBEG))
ENDIF
ENDDO
ENDIF
ENDDO
RETURN
END
FUNCTION POLYL(NL,X)
C NL
is the number of data points in the interpolation
C X
is the coordinate of the output value
C XDAT(NBEG) is the data point just below x
C XDAT(NSKIP) is not included in the interpolation
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NDAT=55,NLAG=12)
DIMENSION ALAG(NLAG),DLAG(NLAG),DDLAG(NLAG)
DIMENSION XDAT(NDAT),FDAT(NDAT)
SAVE NC,XDAT,FDAT
DATA NC/0/
IF(NC.EQ.0)THEN
OPEN(1,FILE='H2.TXT')
DO I=1,NDAT
READ(1,*)XDAT(I),FDAT(I)
ENDDO
CLOSE(1)
ENDIF
IF(NL.GT.NLAG)THEN
PRINT*,' EXCEEDED ALAG DIMENSION IN
POLYL'
READ(*,*)ITEST
STOP
ENDIF
CALL LOCATE(X,XDAT,NDAT,NBEG)
MBEG=MAX0(0,NBEG-NL/2)
MBEG=MIN0(NDAT-NL,MBEG)
NSKIP=0
CALL UELAG(NL,X,MBEG,NSKIP,ALAG,XDAT)
POLYL=0
DO I=1,NL
POLYL=POLYL+ALAG(I)*FDAT(I+MBEG)
ENDDO
RETURN
END
Consider again the BLI integrals between 0 and 1
Use the BLI to find
the best set of points, xdat, and fdat. Then use Gauss Quadrature to find the function
Gdat
Continue with Locate
and Poly to write a function such that
using
Lagrange interpolation.
What is the order of
the error? How would you recommend
estimating the error?