The Lagrange Polynomial (Press Chapter 3)
If we have N distinct data points (xi,fi) , and if we require that the polynomial approximatilon pn(x) be equal to f(x) at all N points. The N constants in the expansion
are uniquely determined by the N equations
which could be written as a set of N linear equations for the c’s. Noting that this implies that there is one and only one polynomial of degree N-1 which passes through the N data points, we can write the Lagrange polynomial which accomplishes this almost by inspection. The Lagrange polynomials are defined with respect to x and the data abscissa xk and xj as
.
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 Lk is a polynomial of degree N-1 and that sums of the Lk 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.
Thus the polynomial
which depends on x only through the Lk(x) terms is explicitly equal to f(xk) for x=xk and of order N-1. When a function is known at a few points over a
region as shown. 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.
4. Coding for the Lagrange Polynomial
A quotation from our text “It is not terribly wrong to implement the Lagrange formula straightforwardly, but it is not terribly right either”[1] gives fair warning that a few extra computer steps are involved in coding the above in directly and that many people will say that it lacks elegance. It is very easy and also easy to debug. In fortran
--- The Fortran code for this is in for\TLagrang.for, this code includes for\lagrange.for and also for\LOCATE.FOR which is described in locate.doc. A much more complex set of code is for\FINDFUN.FOR which is described in Findfun.doc.---
C$NOWAR
C *** FIRST EVALUATE THE FUNCTION AT N DATA POINTS
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NDAT=100)
COMMON /F1/XDAT(NDAT),FDAT(NDAT)
PRINT*,' F(.1),FA(.1),DIFF'
DO N=2,99
BEG=0
H=.2D0
DO I=1,N
XDAT(I)=BEG+H*(I-1)
FDAT(I)=FUN(XDAT(I))
ENDDO
NL=N
X=.1
WRITE(1,*)N,ABS(FUN(X)-POLYL(NL,X))
C5 PRINT*,' ENTER X'
C READ(*,*)X
C IF(X.NE.-1.)THEN
C FA=POLYL(NL,X)
C F=FUN(X)
C PRINT*,FA,F,F-FA
C GOTO 5
C ENDIF
ENDDO
END
FUNCTION POLYL(NL,X)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NDAT=100)
COMMON /F1/XDAT(NDAT),FDAT(NDAT)
DIMENSION ALAG(NDAT)
CALL UELAG(NL,X,ALAG,XDAT)
POLYL=0
DO I=1,NL
POLYL=POLYL+ALAG(I)*FDAT(I)
ENDDO
RETURN
END
SUBROUTINE UELAG(NL,X,ALAG,XDAT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ALAG(NL),XDAT(NL)
DO M=1,NL
ALAG(M)=1.D0
DO J=1,NL
IF(J.NE.M)THEN
ALAG(M)=ALAG(M)*(X-XDAT(J))/(XDAT(M)-XDAT(J))
ENDIF
ENDDO
ENDDO
RETURN
END
FUNCTION FUN(X)
IMPLICIT REAL*8 (A-H,O-Z)
FUN=1/(.5D0+X*X)
RETURN
END
Figure 1 Error in lagrange interpolation as a function of the order of the interpolating polynomial
The test code for the while statement and the formats in c is
--- The C code for this is in cpp\TLagrang.c, this code includes cpp\Lagrange.c and also Locate.c which is described in locate.doc. A much more complex set of code is cpp\findfun.c which is described in Findfun.doc.---
#include <stdio.h>
/* while test */
main()
{
double x;
x=0;
while(x!=-1){
printf("enter a value for x\n");
scanf("%lf",&x); /* ell f is for double*/
printf("you entered %5.3lf\n",x);}
}
The test code for opening files in c is
#include <stdio.h>
#include <stdlib.h>
/* file test */
void main(void)
{
FILE *fp; /* defines the variable fp as a pointer to a file name */
int i;
/* open myfile for output */
if((fp=fopen("myfile","w"))==NULL){printf(" cannot open file\n");
exit(1);}
/* write formatted output to the file */
for(i=0;i<=10;i++) {
fprintf(fp," testing %d\n",i);
printf(" testing %d\n",i);}
fclose(fp); /* close the file */
}
The code lagi.c is
#include <stdio.h>
#include <math.h>
#define H 0.2 /* step size */
struct {double xdat[100], fdat[100];} f1; /* similar to Fortran named common */
double fun(double x); /* prototype function defs, K&R claim helps
to debug functions defined below */
double polyl(int nl,double x); /* functions copy non arrayed arguments, but arrays refer
back to the calling code */
int uelag_(int nl,double x,double alag[100],double xdat[100]);
void main(void)
{
/* first define the elements in the structure */
FILE *fp; /* makes fp a variable that contians the address of a file*/
int n,i,iend;
double beg,x,f,fa,diff;
for (n = 2; n <= 5; ++n) {
beg = 0.;
iend = n;
for (i = 0; i <= iend; ++i) {
f1.xdat[i] = beg + H * (i);
f1.fdat[i] = fun(f1.xdat[i]);
/* printf(" %f %f\n" ,f1.xdat[i],f1.fdat[i]);*/
}
/* } */ FIRST DEBUGGING BREAK POINT
/* then compute the function and its approx at x */
/* while(x!=-1){ USED WHILE TESTING LAG
printf("enter a value for x\n");
scanf("%lf",&x);
f=fun(x);*/
/* printf(" after fun, before polyl\n");*/
/* fa=polyl(n,x);
printf("x, %5.3lf\n",x);
printf(" fa= %15.12lf\n",fa);
printf(" f= %15.12lf\n",f);
printf(" dif %15.12lf\n",fa-f);
}
x=0; */
/* or write the abs diff to a file */
/* open myfile for output */
x=.1;
f=fun(x);
fa=polyl(n,x);
diff=fabs(f-fa);
if(n==2){if((fp=fopen("file1","w"))==NULL)
{printf(" cannot open file\n");
break;}}
fprintf(fp," %d %10.4le\n",n,diff);
printf("f-fa %10.4e\n",diff);
}
fclose(fp);
}
double fun(double x)
{
double ret_val;
ret_val = 1 / (x * x + .5);
return ret_val;
}
double polyl(int nl,double x)
{
int i;
double ret_val;
double alag[100];
/* printf(" in poly before uelag \n");*/
uelag_(nl, x, alag, f1.xdat);
/* printf(" testing alag \n");*/
/* for (i=0;i<=nl;++i){printf(" %lf \n",alag[i]);}*/
ret_val = 0.;
for (i = 0; i <= nl; ++i) {
ret_val += alag[i] * f1.fdat[i];
}
return ret_val;
}
int uelag_(int nl,double x,double alag[100],double xdat[100])
{
int m,i,j;
i=0;
/* printf(" in uelag nl= %d i= %d \n",nl,i);*/
for (m = 0; m <= nl; ++m) {
/* printf(" m= %d ",m);*/
alag[m] = 1.;
/*printf(" after setting alag[m] \n");*/
for (j = 0; j <= nl; ++j) {
if (j != m) {
alag[m] = alag[m] * (x - xdat[j]) / (xdat[m] - xdat[j]);
/*printf(" bottom loop m,j %d %d \n",m,j);*/
} /* after if loop*/
} /* after j loop*/
} /* after m loop*/
/* for (i=0;i<=nl;++i){printf
("in uelag alag %d %lf \n",i,alag[i]);}*/
return 0;
}
5. Assignment.
Write a Lagrangian interpolation code using the closest four points for evenly spaced data and test it on 1/(1+((x-x0)/w)2 (a Lorentzian) .