Bli 
cpp\cloren.zip for\TLORENT.ZIP cpp\CTLAG.ZIP for\LAGRANGE.ZIP
#include <stdio.h>
#include <math.h>
double aloren(double t);
int bli(double xi[],double fi[],double err[],int np,
int ndo, double (*fext)(double));
void main(void)
{double xi[5001],fi[5001],err[5001];
int np=5000,ndo=2,i,ibli,itest;
FILE *fp;
fp=fopen("aloren.out","w");
if(fp ==NULL){printf(" cannot open bli.out file\n");
scanf(" %d ",&itest);
return;}
printf(" enter the beginning and ending x values \n");
scanf(" %lg %lg",&xi[0],&xi[1]);
printf(" xi[0] %lg xi[1] %lg \n",xi[0],xi[1]);
ibli=bli(xi,fi,err,np,ndo,aloren);
for(i=0;i<ibli;i++)fprintf(fp,"%14.6lg %14.6lg\n",xi[i],fi[i]);
fclose(fp);
return;}
double aloren(double xt)
{double x0=3.47,w=0.57,arg,retval;
arg=(xt-x0)/w;
retval=1/(1+arg*arg);
return retval;}
Figure 1 Loren.out
Lagrange Disclaimer
It is so easy to use more and more Lagrange Polynomials to approximate a function that the limits of its accuracy need to be stated at the outset.
w=z-zm
(1.1)
where ξ is some value in the interval between z and zm. Typically the high derivatives of a function eventually begin to rise faster than (N+1)! . Thus eventually the coefficient of w(N+1) becomes large so that the series is asymptotic in its convergence. (at first gets better as N becomes larger, than gets worse). In fairness, it should be stated that the number of terms where this begins is large and that the numbers sort of oscillate about the bottom.
Figure 2 Demonstration of asymptotic convergence FUN=1/(.5D0+X*X)
The Lagrange Polynomial (Press Chapter 3)
If we know a function at N distinct data points and if we require that the polynomial approximation pn(x) be equal to f(x) at all N points. The N constants in the expansion
(1.2)
are uniquely determined by the N equations
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
One term is always left out of the product that covers N values in which x always appears once. Thus 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 numerator 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 desired is
Note
This depends on x only through the Lk(x) terms. At x = xk, the L’s are equal to zero except for the one multiplying fk which is equal to 1. Thus pn(x) goes through all of the data points and is explicitly of order N-1.
Equation (1.5) is the solution to the problem leading to (1.3). Note that (1.5) does not explicitly give us the coefficients c0 through cN-1. This is a drawback.
The coding for (1.6) and higher derivatives along with a practical example is in ..\wsteve\FDERS\The Lagrange Polynomial.htm
Coding 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”[2] 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 for/lagrange.for -- for\LAGRANGE.ZIP
FUNCTION POLYL(NL,NDAT,X,XDAT,FDAT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XDAT(*),FDAT(*)
DIMENSION ALAG(24)
IF(NL.GT.24)THEN
PRINT*,' FROM POLYL ALAG DIMENSION EXCEEDED NL = ',NL
NL=24
READ(*,*)ITEST
ENDIF
CALL LOCATE(X,XDAT,NDAT,JB)
C*** JB is such that XP(JB)<X<XP(JB+1)
C*** MOVE BACK (NL-2)/2 IF POSSIBLE
C*** MOVE BACK TO LEAVE NL+JB=NDAT
JB=MIN(JB,NDAT-NL)
JB=MAX(0,JB-1-(NL-2)/2)
CALL UELAG(NL,JB,X,ALAG,XDAT)
POLYL=0
DO I=1,NL
POLYL=POLYL+ALAG(I)*FDAT(I+JB)
ENDDO
RETURN
END
SUBROUTINE UELAG(NL,JB,X,ALAG,XDAT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ALAG(NL),XDAT(*)
DO M=1,NL
ALAG(M)=1.D0
DO J=1,NL
IF(J.NE.M)THEN
FACT=(X-XDAT(J+JB))/(XDAT(M+JB)-XDAT(J+JB))
ALAG(M)=ALAG(M)*FACT
ENDIF
ENDDO
ENDDO
RETURN
END
cpp\CTLAG.ZIP
double laglo(double x)
{static int nc=0, ndat;
static double xp[137],dp[137];
FILE *fp;
if(nc == 0)
{nc=1;
fp = fopen("loren.out","r");
if(fp == NULL)
{printf(" cannot open file loren.out \n");
return 0.0;}
ndat=0;
redbeg:
if(fscanf(fp, "%lg %lg", &xp[ndat],&dp[ndat])!= EOF)
{/*printf("The line read was: %lg %g\n", xp[ndat],dp[ndat]);*/
ndat = ndat+1;
if(ndat <= 137) goto redbeg;
else
{printf(" from laglo ndat > 137, enter to continue ");
scanf(" \n");}}
fclose(fp);
ndat=ndat-1;}
return bpolyl(4,ndat,x,xp,dp);}
void main(void)
{double xi[502],fi[502],err[502];
int np=500,ndo=2,i,ibli,itest;
/* double x,test;
double x=1.99,test; */
FILE *fp;
/* test=laglo(x);
printf(" x, TEST =%lg %lg \n",x,test);
return; */
fp=fopen("blic.out","w");
if(fp ==NULL){printf(" cannot open bli.out file\n");
scanf(" %d ",&itest);
return;}
printf(" enter the beginning and ending x values \n");
scanf(" %lg %lg",&xi[0],&xi[1]);
printf(" xi[0] %lg xi[1] %lg \n",xi[0],xi[1]);
/* for(i=0;i<=np;i++)
{x=-4+i*(xi[1]-xi[0])/np;
test=laglo(x);
fprintf(fp,"%14.6lg %14.6lg\n",x,test);} */
ibli=bli(xi,fi,err,np,ndo,laglo);
/* printf(" after bli xi[0] = %lg \n",xi[0]); */
for(i=0;i<=ibli;i++)fprintf(fp,"%14.6lg %14.6lg\n",xi[i],fi[i]);
fclose(fp);
return;}
#include <stdio.h>
int locate(double x,double r[],int n);
void uelag(int nl,int jb,double x,double alag[],double xdat[])
{int m,j;
double fact;
for (m = 0; m < nl; ++m)
{alag[m] = 1.;
for (j = 0; j < nl; ++j)
{if (j != m)
{fact = (x - xdat[j+jb]) / (xdat[m+jb] - xdat[j+jb]);
/* printf(" from uelag fact = %lg \n",fact); */
alag[m] = alag[m] * fact;}}}
return;}
double bpolyl(int nl,int ndat,double x,double xdat[],double fdat[])
{int i,jb;
double ret_val;
double alag[24];
if(nl > 24)
{printf(" from polyl alag dimension exceeded nl = %d \n",nl);
nl = 24;
scanf(" \n");}
jb = locate(x,xdat,ndat);
if(jb > ndat - nl)jb = ndat - nl;
jb=jb-(nl-2)/2;
if(jb < -1 ) jb = -1;
uelag(nl,jb,x,alag,xdat);
ret_val = 0.;
for (i = 0; i < nl; ++i)
{ret_val += alag[i] * fdat[i+jb];}
return ret_val;}
The one change from defaults is to switch from full debugging information to line number information.
#include <stdio.h>
#include <math.h>
These are needed for all files in Watcom, not for most other compilers.
double bpolyl(int nl,int ndat,double x,double xdat[],double fdat[]);
int bli(double xi[],double fi[],double err[],int np,
int ndo, double (*fext)(double));
/* #include "d:\cpp\bli.c"
#include "d:\cpp\lagrange.c"
#include "d:\cpp\locate.c" */
These are commented out here, they are needed for other compilers.