Assign2

Bli Integration 2 The Lagrange Polynomial

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

Assignment 2

            Consider again the BLI integrals between 0 and 1

  1. ftest = x4
  2.  Integral should be 0.886226925455
  3.  Integral 0 to 1 should be 1.64493406685
  4.   error integral.

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?