Best Linear Interpolator for\bli.zip cpp\cbli.zip
The routine starts with the 4 red points. Points 1 and 3 are interpolated to find a value for point 2. Points 2 and 4 are interpolated to find a value for point 3. The region with the largest value of d times delta x has a point placed in the middle of it. This continues until the maximum number of points is reached or until the maximum interpolation error falls below some tolerance.
The object is to use exactly N function evaluations to get the best estimate of a function possible with N function calls. - The assumption is that each call is relatively time consuming.
Figure 1 Region of a curve showing linear interpolations for points 3 and 6.
The BLI code
The code begins with a small number of points as shown in figure 1. Each point is assigned a relative error based on the difference between its value and a linear interpolation between the neighboring points. The error at the beginning and ending points are taken to be those at the second and second to last points. Each interval is assigned an error based on the absolute value of the sum of the end points of the interval times the size of the interval to a power. The interval with the largest error has a value of the underlying function evaluated at its mid point. In figure 1, this would involve placing a new point in the interval from 2 to 3. The new point changes the error estimates for points 1, 2, new, and 3. All others remain the same.
Calling code
Fortran for\TBLI.FOR
IMPLICIT REAL*8 (A-H,O-Z) single precision is almost never adequate for
physics
DIMENSION XI(5001),YI(5001),ERR(5001) gplot can handle 5001 points
EXTERNAL ALJ tells
the compiler that ALJ is a function rather than an array
NP=5001
PRINT*,' ENTER B,E'
READ(*,*)XI(1),XI(2)
NDO=2 This tells the BLI subroutine to use
the points in a uniform grid, if ndo ≠ 2. These
are dimensioned variables because the colling routine frequently knows about
peaks valleys and special points. These
can be put in XI(I) making NDO > 2.
When this happens the BLI will not begin by putting half of its points
in a uniform mesh.
CALL BLI(XI,YI,ERR,NP,NDO,ALJ)
OPEN(1,FILE='BLI.OUT')
WRITE(1,'(2G14.6)')(XI(I),YI(I),I=1,NP) the format is single precision for gplot.
STOP
END
C cpp\tbli.c
#include <stdio.h>
double alj(double t); function prototype
int bli(double xi[],double fi[],double err[],
int np,int ndo,double
(*fext)(double)); function prototype continues over two lines
note
the signal that the function is passed.
void main(void)
{double xi[5001],fi[5001],err[5001];
int
np=511,ndo=2,i,ibli,itest; code uses 512 points from 0 to 511, ndo ≠ 2, bli
expects to find ndo initial values of xi and fi.
FILE *fp; allows fp to contain file information
fp=fopen("blic.out","w"); attempts to open output file
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]); note the &’s
scanf is a cross that C programmers bear. The & causes the result from scanf to be
placed in the addresses of the array elements.
There is always some doubt as to how array elements are passed. Removing the &’s produces a satisfying
core dump.
printf(" xi[0] %lg xi[1] %lg \n",xi[0],xi[1]);
ibli=bli(xi,fi,err,np,ndo,alj); since bli is a function, it returns the number
of points, which could also have been np
for(i=0;i<ibli;i++)fprintf(fp,"%14.6lg
%14.6lg\n",xi[i],fi[i]); note the formatted 14.6 to output single
precision. The 6 in the format comes about because it is the effective limit for
the single precision in gplot. The 14 is
because with all of the overhead of formatted output, this is what it takes.
See Format.txt
fclose(fp);
return;}
Beginning of BLI
Fortran for\BLI.FOR
SUBROUTINE BLI(XI,FI,ERR,NP,NDO,FEXT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ERR(*),XI(*),FI(*)
DATA ALPHA,BETA/1D-1,1D-2/ see #Error
determination
C *** USE HALF THE POINTS IN AN EVENLY SPACED MESH
IF(NDO.EQ.2)THEN NDO in this routine will be the current number
of points
NDO=NP/2
H=(XI(2)-XI(1))/(NDO-1)
DO I=1,NDO
XI(I)=XI(1)+(I-1)*H
FI(I)=FEXT(XI(I))
ENDDO
ENDIF
C
#include <math.h>
int bli(double xi[],double fi[],double err[],int np,
int ndo, double (*fext)(double))
{int i,ib,ie,itest,ile,iadd;
double h,errl,errc,errm,xb,xe,fint;
double alpha=0.1,beta=0.01,div; see #Error
determination
if(ndo == 2)
/* *** USE HALF THE POINTS IN AN EVENLY SPACED MESH
note that C starts with the point at 0, but still has ndo points and
ndo-1 intervals */
{ndo=np/2;
h=(xi[1]-xi[0])/(ndo-1);
for (i=0;i < ndo;++i)
{xi[i]=xi[0]+h*i;
fi[i]=(*fext)(xi[i]);}}
Error determination
The 8th power of the interval used
in previous years is being replaced by a limited error size.
When functions vary over many orders of magnitude, the error on a point is not simply the distance from the interpolated value to the actual point. The relative error can be the more relevant one. In curve fit applications (../Const/denton/stanfit/Stanfit3.html#Practical) it was found appropriate to use an expression similar to
The minimization is then of
In BLI the appropriate error at each point is
The term β is naturally set to 0.01 for graphical applications, while the tendency of physics functions to be on the oreder of 1, implies that alpha should initially be 0.1.
Fortran 
C *** NOTE THAT THERE ARE ONLY NDO-1 INTERVALS
C *** ERR(I) REFERS TO THE INTERVAL FROM X(I) TO X(I+1)
DO I=2,NDO-1
XB=XI(I-1)
XE=XI(I+1)
IF(XB.EQ.XE)THEN
FINT=.5D0*(FI(I-1)+FI(I+1))
ELSE
FINT=FI(I-1)+((XI(I)-XB)/(XE-XB))*(FI(I+1)-FI(I-1))
ENDIF
DIV=BETA*ABS(FI(I))
IF(DIV.LT.ALPHA)DIV=ALPHA
ERRC=ABS(FINT-FI(I))/DIV
IF(I.EQ.2)ERRL=ERRC
ERR(I-1)=(XI(I)-XI(I-1))*(ERRC+ERRL)
ERRL=ERRC
ENDDO
ERR(NDO-1)=2*(XI(NDO)-XI(NDO-1))*ERRL
C
for (i=1;i < ndo-1;++i)
there are points up to ndo-1
The loop usesfi( i+1) implying that this loop
canonly go to ndo-1
{xb=xi[i-1];
xe=xi[i+1];
fint=fi[i-1]+((xi[i]-xb)/(xe-xb))*fi[i+1];
div=beta*fabs(fi[i]);
if(div<alpha)div=alpha;
errc=fabs(fint-fi[i])/div;
if(i==1)errl=errc; needed for first interval
err[i-1]=(xi[i]-xi[i-1])*(errc+errl);
This will end with i=ndo-2 and thus find all
but the last interval
errl=errc;}
err[ndo-2]=(xi[ndo-1]-xi[ndo-2])*errl; there are points up to ndo-1, but intervals
only to ndo-2
Fortran
50 CONTINUE
This address will be used after each point has
been added
ERRM=ERR(1)
ILE=1
DO I=2,NDO-1 < There are intervals at 1,2, …, NDO-1 for NDO points
IF(ERR(I).GT.ERRM)THEN
ERRM=ERR(I)
ILE=I
ENDIF
ENDDO
C
f_largest_err: This is a label, greatly discouraged in modern
programming, it allows this segment to respond as each ndo point is added. Assignment
make this into a do loop and test resulting
code.
errm=err[0];
ile=0;
for
(i=1;i<ndo-1;++i) There are points at 0,1,2,…,ndo-2
{if(err[i]>errm){errm=err[i];
ile=i;}}
/* printf(" largest err is %d %lg \n",ile,errm);
scanf("
%d",&itest);*/ compiler will usually comment on itest, it is
used here in testing code.
Fortran
NDO is not incremented until the end of this section. Thus for most of the section there are NDO+1 points.
C *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT ILE+1 BECOMES
C *** ILE+2
XI(NDO+1)=XI(NDO)
FI(NDO+1)=FI(NDO) these are added outside the do loop because
there is no err(ndo)
DO I=NDO-1,ILE+1,-1 note that the do loop is decremented to keep
from overwriting terms
XI(I+1)=XI(I)
FI(I+1)=FI(I)
ERR(I+1)=ERR(I)
ENDDO
IADD=ILE+1
XI(IADD)=(XI(ILE)+XI(ILE+1))/2
FI(IADD)=FEXT(XI(IADD))
C
/* *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT ILE+1 BECOMES ILE+2 */
xi[ndo]=xi[ndo-1]
fi[ndo]=fi[ndo-1]
for (i=ndo-2;i>=ile+1;--i)
{xi[i+1]=xi[i];
fi[i+1]=fi[i];
err[i+1]=err[i];}
iadd=ile+1;
xi[iadd]=(xi[ile]+xi[ile+1])/2;
fi[iadd]=(*fext)(xi[iadd]);
Recalculating error estimates
Figure 2 A new point inserted into interval 2
In figure 2, it is assumed that the the old interval 2 was found to have the largest error estimate. A new point is added at the x. This results in different fint values for points 2, x, and 3. This changes the error estimate for intervals 1, 2, 3 and 4, but intervals beyond 5 are unchanged.
Fortran
C *** ERRL FOR ILE-1 AND ERR(ILE-1 TO ILE+1) NEEDS TO BE RECALCULATED
IB=MAX(2,ILE-2)
IE=MIN(NDO,ILE+3)
DO I=IB,IE
XB=XI(I-1)
XE=XI(I+1)
FINT=.5D0*(FI(I+1)+FI(I-1))
IF(XE.GT.XB)THEN
FINT=FI(I-1)+((XI(I)-XB)/(XE-XB))*(FI(I+1)-FI(I-1))
ENDIF
DIV=BETA*ABS(FI(I))
IF(DIV.LT.ALPHA)DIV=ALPHA
ERRC=ABS(FINT-FI(I))/DIV
IF(I.EQ.2)ERRL=ERRC
IF(I.GT.ILE-1)ERR(I-1)=(XI(I)-XI(I-1))*(ERRC+ERRL)
ERRL=ERRC
ENDDO
IF(ILE+3.GE.NDO)THEN
ERR(NDO)=2*(XI(NDO)-XI(NDO-1))*ERRL
ENDIF
NDO=NDO+1
IF(NP.GT.NDO)GOTO 50
RETURN
END
C
/* *** ERRL FOR ILE-1 AND ERR(ILE-1 TO ILE+1) NEED TO BE RECALCULATED */
ib=ile-1;
if(ib<1)ib=1;
ie=ile+3;
if(ie>ndo-1)ie=ndo-1;
for(i=ib;i<=ie;++i)
{xb=xi[i-1];
xe=xi[i+1];
fint=fi[i-1]+((xi[i]-xb)/(xe-xb))*fi[i+1];
div=beta*fabs(fi[i]);
if(div<alpha)div=alpha;
errc=fabs(fint-fi[i])/div;
if(i==1)errl=errc;
if(i>ile-1)err[i-1]=(xi[i]-xi[i-1]*(errc+errl);
errl=errc;}
if(ile+3>=ndo-1)err[ndo-1]=2*(xi[ndo-1]-xi[ndo-2])*errl;
ndo=ndo+1;
if(np > ndo)goto f_largest_err;
return ndo;}
Figure 3 Points for Lennard Jones potential between 0 and 100 256 points
Figure 4 Detail on linear scale of interesting region of Lennard Jones potential using 256 points.
scanf(" %lg %lg",&xi[1],&xi[2]);
printf(" xi[1] %lg xi[2] %lg \n",xi[1],xi[2]);
ibli=bli(xi,fi,err,np,ndo,alj); Since bli is a function, it returns a
value, normally ibli is equal to np, but for debugging, one may want to stop
early. Note that the name alj is passed
as simply another argument. The fact
that it is a function is above.
for(i=0;i<=ibli;i++)fprintf(fp,"%
14.6lg %14.6lg\n",xi[i],fi[i]); note the formatting (See Format.txt) and the \n for linefeed
fclose(fp);
return;}
Consider the Lennard Jones Potential .doc
The
anti-symmetry of parallel spin electrons makes them truly resent being pushed
together.
The slight polarization of the electrons in one atom relative to another gives rise to a slight attraction between the atoms.
The Lennard Jones potential
Fortran for\LENNJONE.FOR
FUNCTION ALJ(XT)
IMPLICIT REAL*8 (A-H,O-Z)
X=MAX(1D-2,XT) one always has to watch out for
overflows. Do not change xt
the change will also be in the main routine.
ALJ=1/X**12-1/X**6
RETURN
END
C cpp\lennjone.c
double alj(double xt)
{double xi,xi3,xi6,xi12,retval;
xi=xt; This is not Fortran xt is
a copy of the variable in the main code that is not returned. - I did not need to use xi.
if(xt<.01)xi=.01;
xi3=1/(xi*xi*xi);
xi6=xi3*xi3;
xi12=xi6*xi6;
retval=xi12-xi6;
return retval;}